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BY 
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1. Introduction. It has been observed in the irrigation wells of the Hawaiian Islands 
that the water-level fluctuations have frequency components corresponding to those 
of the ocean tides [1]. This phenomenon was analyzed by Carrier and Munk [2], assuming 
the observed ground-water fluctuations to represent a diffusive transmission of the 
tidal disturbances through the porous voleanic structure of the island. The purpose of 
the investigation was to use the results in estimating the permeability of the porous 
medium. 

In [2] it was assumed that the porous medium was infinitely deep. In actual fact, how- 
ever, there will be an essentially impenetrable bounding surface (see Fig. 1). This paper is 
concerned with the analysis of the same problem treated in [2], but taking account of 
the bounding bottom surface. Numerical computations are carried out for several values 
of the dimensionless depth. Also the limiting case of shallow water theory is studied. 
Using the results of the infinite depth, shallow depth, and finite depth theory it is possible 
from the graphs given in Figs. 2 and 3 to estimate the amplitude and phase lag in the 
fluctuations of the ground-water as a function of the distance inland for various values 
of the dimensionless depth. It is found that for the values of the physical parameters 
which are probably of most concern the infinite depth theory gives satisfactory results 
in the region of interest. 

2. Formulation of the problem. Although the formulation and the first part of the 
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Fig. 2. Ground-water amplitude vs. distance inland, x. (Dotted lines indicate the expected correction 
to the computed curves when all modes are considered) 


analysis of this problem follows quite closely that related in [2] it is convenient, for the 
sake of completeness, to repeat part of that work here. As the water in the ocean bounded 
by Oy’ and OB (see Fig. 1) rises and falls about its mean level OA, the pressure on the 
line OB varies. Corresponding to this periodic change in pressure on OB we can expect 
periodic fluctuations in the free surface of the ground-water, i.e. Ox’. 

The equations governing the motion of the fluid in the porous medium are the con- 


servation of mass 
> fa) io 
div (pv) = 5p (PY) (2.1) 
( 


and Darcy’s law which replaces the conservation of momentum law (see [3]), 

v = —(k/u) grad (p — po). (2.2) 
Here p and uz are the density and the viscosity of the fluid; @ and k are the porosity and 
the permeability of the material; v is the velocity with components u and v in the 2’ 
and »’ direction; p is the gauge pressure; pp = —pogy’ Is the pressure when the fluid 
motion is zero; pp is the mean density; and subscript notation indicates partial differ- 
entiation. The simple compressibility law used by Carrier and Munk is 


pO = poAi1 + Sp — po}; (2.3) 


where 4 is essentially (p c”)~* with c the speed of sound in the fluid. It should be pointed 
out that Eq. (2.2) says that the pressure gradient is proportional to a velocity rather 
than an acceleration as in the Navier-Stokes equation. As a result we will obtain finally 
an equation of the diffusion type rather than a wave equation; hence the free surface 


; 


amplitude will decay in x 
rr . ° ’ 
The boundary conditions expressed in terms of the pressure p are 
p=0, onthe free surface, (2.4a) 


y= — poJYon Tr Ue on OB, (2.4b) 
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a(p — p.)/dy’ = Oony’ = —H, (2.4c) 


p = 0, on the free surface, 

where H is the depth of the porous medium. The last boundary condition results from 

the requirement that the normal component of velocity be zero on the impenetrable 

bottom. The pressure q, is of course directly proportional to the tidal-wave amplitude. 
If we let g = p — po, and combine Eqs. (2.1), (2.2) and (2.3) we obtain 


Aq = (u68/k)a , (2.5) 


where A is the Laplacian operator. If we denote by n(z’, ¢) the y’ coordinate of the free 
surface, Eq. (2.4a) implies q(x’, n, 1) = pogn, but on the free surface 4, = v/e, hence 
using Eq. (2.2) our boundary condition (2.4a) may be expressed as 


de + (pogk/u®)qy =O on y’ = 7, x>O0. (2.6a) 

Actually as in the usual linear theory of water waves this boundary condition is to be 
applied on y = 0. The boundary conditions (2.4b) and (2.4c) may be written as 

q(on OB) = qe***, (2.6b) 

q, = 0,y’ = -H. (2.6c) 


We shall only solve this problem in the case that the line OB occupies the half-line 
y = 0, x < 0. That is we take ¢ = 0° (see Fig. 1). Actually this is fairly realistic since 
¢ is probably of the order of 5° or so. 

Finally if we introduce the following dimensionless variables 


r=, z= 2'/L, y=y'/L, h = H/L, 


9 
L = (pogk)/(uOow) , € = (pog7kd)/(uOw), q = nel, ye’, sie 
we obtain 

Ag — iep = O, (2.8) 

with the boundary conditions 
¢, tig = 0, y = 0, z> 0, (2.9a) 
g= 1, y = 0, z <0, (2.9b) 
y, = 0, y = —h, —-o <4< om, (2.9¢e) 


The free surface n(x, ¢) is (pog)~’ ¢(2, 0, #), but from Eqs. (2.7) ¢ = q.¢(2, y) exp (twt) so 


n(x, t) = ~ g(x, O)e'*' x > 0. (2.10) 
So the problem of determining the free surface is exactly that of determining ¢(z, 0). 
The combination of parameters q,/pog is the maximum height of the tidal-wave measured 
from y’ == §, 

Before proceeding to a solution of the problem defined by Eqs. (2.8) and (2.9) it is 
perhaps worthwhile to mention briefly the size of the parameters which appear in this 
problem. We have p/p) = 0 (10~*cem*/sec), k = 0(5 X 10°-°cm*), 6 = 0(.20), g = 980 
em/sec’, and w for a twenty-four hour tide is 27/24 hours, hence L = 0(1000 ft). Since 
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ce is 0(5000 ft/sec) for water, « is a very small number, 0(10~*). Finally, a reasonable 
value for the depth of the ocean is about three miles so h may be as large as 15. 

3. Shallow water theory. Before considering the general problem given by Eqs. 
(2.8) and (2.9) let us look at the limiting case in which the depth H is small enough that 
we can neglect variations in the y’ direction, and also set v = 0. Then u(z’, ¢) represents 
an averaged velocity across the section —H < y’ < 0. If we assume incompressibility, 
i.e. 6 = 0, the conservation of mass equation appropriate to this situation is 


Hpwtzs = —poPon .« (3.1) 

Darcy’s equation reduces to 
u= —(k/u)q . (3.2) 
Since there is no variation in the y direction our condition that g = p gy on the free 


surface must hold throughout the strip —H < y’ < 0,2 > 0. Using this and Eqs. (3.1) 
and (3.2) we obtain 
Qerz' — (uO./kpogH)q. = 0, o > 6. (3.3) 


The condition that gq = q, exp (iwt) for x < 0 is now applied at t = 0; consequently 
we set g = 9,¢(x’) exp (wt). Then Eq. (3.3) becomes 


¢'2' — (HL) 'e = 0, 2” > ©. (3.4) 
An appropriate solution of Eq. (3.4) satisfying a finiteness condition at infinity is 


g(x’, t) = exp [— (i/HL)'z’]. 


Hence 
n(x’, t) = (q:/pog) exp {— 2’ (2HL)* + ifwt — a’ /(2HL)*}}, (3 .5a) 
= (41 Pod) exp i— Z (2h)? + t[wt = (2h)*]}. (3.5b) 


Actually, in order for this theory to be valid not only must the wave length of the 
disturbance be large compared to H as in the usual shallow water theory but also H 
must be small compared to the other natural length scale, L, which appears in the 
problem, i.e. h must be small. This can be seen by an examination of the behaviour of 
the solution of the general problem. This is done in Sec. 5 where it is found that for 
h < 1/4 we can expect the shallow water theory to be quite accurate. The amplitude 
and phase lag of pogn(z, t)/q, are plotted in Figs. 2 and 3 as a function of z for h = 1/4. 

4. Solution of the problem. To solve the problem defined by Eqs. (2.8) and (2.9) 
we shall use the method of Fourier transforms and the Wiener-Hopf technique. Let 


(E, y) = | e “g(x, y) dz. (4.1) 


Then the transform of Eq. (2.8) is 


&, — (° + id = 0. (4.2) 


A solution of this equation satisfying the boundary condition (2.9c) is 


GE, y) = AE) cosh {(y + h)C}, (4.3) 


where C = (é’ + ie)'”’ and A(E) is to be determined by satisfying the remaining boundary 


conditions. 
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Let 
where (4.4a) 
Or>0 
tae  *se. (4.4b) 
lyo(x, 0) x > 0 
and 
f(x) = ¢,(x, 0) + i¢(z, 0). (4.5) 
It is clear that o(z, 0) = g,(x) + g2(x); hence* 
P(E, 0) = A(é) cosh Ch = G,() + G.(§), (4.6) 
= (a — it)" + G,@). 
Also using Eq. (4.3) 
F(é) = (C' sinh Ch + 7 cosh Ch) A(é). (4.7) 
Combining Eqs. (4.6) and (4.7) we obtain 
Fé) = K®{G.@ + G.©}, (4.8) 
where 
-,, __ Csinh Ch + 7 cosh Ch 
ko = cosh Ch (4.9) 


If we recall from Sec. 2 that we wish to determine ¢(z, 0) it is clear from Eq. (4.4) 
that our problem is now that of determining g,(z) and hence G,(¢). To determine G,(é) 
using Eq. (4.8) we shall use the Wiener-Hopf technique. This technique has been used 
to treat similar problems (see for instance [2, 4, 5]); consequently the analysis will only 
be briefly outlined here. First G,(é) is analytic in the upper half plane, (UHP), Im(é) 
> — a; G,(£) is analytic in the LHP and F(é) is analytic in the UHP. The function 
K(é) is analytic and non-vanishing in a strip containing the real axis. This will be seen 
clearly at the end of this section where K(é) is represented as the quotient of two infinite 
products. It might be noted that though C = (¢ + 7e)'” is a multivalued function, 
K(&) as defined by Eq. (4.9) consists only of even terms and hence does not have any 
branch points. Assuming for the moment that we can write K(£) as K_(£)/K,(&) where 
K_(£) is analytic and non-vanishing in the LHP, and K,(é) is analytic and non-vanishing 
in the UHP we can rewrite Eq. (4.8) as 

F®@)K.@) — KA(— ia)G,@ = {K-@ — K(— ia)}G,@ + K-©G,(). (4.10) 
The left hand side of this equation is analytic in the UHP, the right hand side is analytic 
in the LHP and they agree in a common strip of analyticity. Hence Eq. (4.10) defines 
an entire function Z(£). We shall show shortly that K_(€) = 0(¢'””) as| £| ~ ©, Im(é) < 0 
and K,(£) = 0(¢°'”*) as| £ | — «, Im(é) > 0. Using this and investigating order con- 
ditions at infinity we can show that £(é) = 0, consequently 


G,(f) = (Kt io -- i} G,(é). (4.11) 


*Capital letters are used to denote the Fourier transform. 
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It is now necessary to determine K_(£) and K,(£). The splitting of K(&) is done in 
a manner exactly analogous to that used by Heins in [4] and [5]. Using the infinite 
product representation of cosh z, (see [6]) we have 


=) 


M(®) = cosh Ch = J] {1 + (2Ch)*/(2n + 1)*x"} = m@)m(— 8, = (4.12) 
where 
m(t) = TI {{1 + (4ieh®)/(Qn + 1)*x?} (4.13) 


+ i2hé/(2n + 1)r} exp [— i2hE/(2n + 1)r]. 


We have inserted the exponentials to insure absolute convergence of the infinite products 
defining m(£) and m(—é) in the UHP and LHP respectively. If we write M(£) as M_(&)/ 
M.(é) it is clear that M_(&) = m(£) has no zeros or poles in the LHP. Similarly 1/4 , () 
m(—é£) has no zeros or poles iy the UHP. 

The function L(é) = C sinh Ch + 7 cosh Ch has zeros at Ch = + 78,,n = 0,1, 2,---, 
where the 8, are complex numbers lying in the first quadrant. For n large they may be 
determined by the asymptotic relation 8, = nr + th/nx + 0*([nz]’). We may write 
L(é) as 


iT] {1 + (Ch)?/B2} = ilU— 8, (4.14) 


Fee 


where 


+ (teh?) /B2)? + ith/Bo} {{1 + (teh?) /B2]? + ith/B,} 
l s f 5 


n=1 


= 
rr 
I 


(4.15) 
exp (— t&h/nm). 


Again we have inserted the exponentials in order to insure absolute convergence in the 
appropriate half planes. If we write L(é) as L_(£)/L,(&) and take L_(&) = U(&), 1/L.(&) = 
il(—£) it is clear that L_() is free of zeros and poles in the LHP and L.(é) is free of 
zeros and poles in the UHP. 

Consequently we have 


K_(é exp {x(&)}L_(6)/M_(é) = exp {x(é) j l(&)/m(&), (4.16a) 
K .(é exp {x(£)}L.(€)/M.(&) = exp {x(&)}m(— )/il(— &). (4.16b) 


We shall choose the factor exp {x(£)} introduced in Eqs. (4.16a) and (4.16b) in such a 
] g | \ 


manner that K =< and K.(¢) have algebraic behaviour as | | — ~ in the LHP and 


UHP respectiv 


To in} vestigate ie behay iour of K (£) for Im(é) < & E — @& We first note that 
the terms involving « may be neglected against unity for | | — ©. Since 6, — nm 


asn— we have that K_(é) is of the order 
exp {x(&)}(1 + 7&h/Bo) I] {(1 + w/n) exp (— w/n)} / 


I] (11 + 2w/(2n + 1)] exp [— 2w/(2n + 1)]} 











n 
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where w = 7th/x. Now using the relation that 


1/T(w) = we™” [J] (1 + w/ne™*”, 


n=1 


and Stirling’s asymptotic formula for the gamma function, (see [6]) we obtain 


K_(¢) = 0{w' exp [x(@) + win4]}, 


for Im(é) < 0,| | — &. So choosing x(¢) = — wln 4 = — (ith ln 4)/m we have that 
K_(£) = 0(€'*) as | & | — , Im(£) < 0. A similar argument will show that K,(é) = 
O(w '’*) for | £| — @, Im(£) > 0. With these order relations it is not difficult to show 


that E(¢) is zero as mentioned earlier; and hence we obtain Eq. (4.11) for G,() where 
K_(£) is defined by Eq. (4.16a). In particular it can be seen from Eq. (4.16a) and the 
definitions of 1(€) and m(£é) that K_(0) = 1. 

Using the usual inversion formula we have that 


) 


G.(x) = (2m)™' | e'*{(K_(— ia) — K_()]/(a — i&)K_(€)} dé, (4.17a) 


ial f ‘ = K_(— ia) cosh Ch — K.(&)K(é) 
= (27 s (a — 2&)K,(E) K(E) 


In the limit as a > 0 it is clear from Eq. (4.17a), that G,(£) is not singular at the origin; 
hence we may actually take the real axis as our path of integration in evaluating g.(z). 
Of course we shall actually close our path of integration in the UHP when z > 0 and 
in the LHP when x < 0. For the case x > 0 it is convenient to use Eq. (4.17b) for evaluat- 
ing g.(x). Since K..(é) is non-vanishing in the UHP, g,(x) will be simply the sum of the 
residues at the poles of the integrand which occur at Ch = 78, ,n = 0,1, --- . Carrying 
out this straightforward computation we obtain in the limit 





dé. (4.17b) 


(4.18) 





~-1 wo (Br "on 
g(x) = h Qe . B2\R (7..\? 
, <= \a,/ (B, + th — h)K,(ia,) 
where a, = (82/h? + ie)'’”. In the particular case that the fluid is incompressible, i.e. 
e = 0, B,/a, = h and Eq. (4.18) becomes 


Qnz 


P 
(rt) = h < x Saree eS “Se = 8, / . =. 
G2\%) é' dX (hea, + th — W)K.(ia,)’ ™ B./h (.-9) 


5. Numerical computations and discussion. In this section we shall only be con- 
cerned with the case in which the fluid may be considered as incompressible, then g,(x) 


is given by Eq. (4.19). In [2] a few values of g.(x) were computed for « = .01 and com- 
pared to the « 0 case; the amplitude and phase lag in the ground-water fluctuation 
for « = .01 were slightly lower than for « = 0. Since ¢ is however 0(10~*) we should 


expect very little error in actually setting « = 0. 
irst let us determine when the shallow water theory solution given by Eqs. (3.5) 
may be expected to be valid. In order for g(x) as given by Eq. (4.19) to agree with the 
shallow water solution it is necessary that a) ~ (z ‘h)'”? as h > 0, and also that all the 
coefficients of the higher order terms' must approach zero. It can be shown with little 
\We shall refer to the term exp (—aor) in Eq. (4.19), which is the dominating term as z > ©, as 
the fundamental term. 








336 R. C. DIPRIMA [Vol. XV, No. 4 
difficulty from an investigation of the transcendental equation C sinh Ch + 7 cosh 
Ch = 0 that for small h, 8, ~ (h)'” and hence ay ~ (i/h)'”*. Also upon noting that for 
small h, 8, ~ nx for n > 1, it can be seen from Eq. (4.19) with the aid of the represen- 
tation of K,(za,) given in Eq. (5.1) that all the coefficients of the higher order terms do 
approach zero as h — 0. Hence g.(x) as given by Eq. (4.19) does approach the shallow 
water solution as h — 0. To determine quantitatively when Eq. (3.5b) is valid we have 
computed 6, and a, as a function of h. Also the ratios of the wave length predicted by 
the shallow water theory, \ = 2r/(2HL)'”’, to H and that of the fundamental mode, 
Xo = 2x/Im(ay), to H have been computed. These results are given in Table 1 and 


TABLE 1. 
h is the dimensionless depth, \ = 2x(2HL)*? is the wave length for shallow water theory, \x = 2xL/Im(ao) 
is the wave length of the fundamental mode for finite bottom theory. 





| | 





h ay(e = 0) \/H | ho/H 
— 0 | — (ih)il2 | 

25 .3676 +7 .3382 1.4704 + 7 1.3528 76 18.57 

.50 .5376 +7 .4548 1.0752 +2 .9096 12.57 13.81 
1.00 £8004 +i .5702 | .8004 +i .5702 8.88 | 11.09 
2.00 1.1828 +7 .5832 5914 +7 .2916 6.28 10.77 
3.00 1.3739 +i .4775 .4579 +4 .1592 185 =| 13.16 
5.00 1.5033 +i .3090 .3007 +7 .0618 3.97 | 20.03 
10.00 1.5547 +7 .1569 1555 +i .0157 2.81 10.04 
15.00 1.5638 +2 .1046 10438 +2 .0070 2.29 59.83 





graphically in Fig. 4. It appears from Fig. 4 that we may expect the shallow water 
theory to be accurate over the entire range of x forh < 1/4. 

In order to compute g.(r) for various values of h it is necessary to cast K, (ta 
a form more suitable for numerical analysis than that given by Eq. (4.16a). This can 
be done in a straightforward manner by using the infinite product representation of the 


into 


gamma function. We obtain, when e = 0, that 


\ 


B,{1(6,/m)]? exp [(8,]n4) /7] , ee a 7 Py 
Drill TA\TM9R /r) (1 + B,/Bm)/( 8 / mr (5.1 
Zu ] + B, /Bo) T'(26,, 7) I] (1 . . \ + f m Ie ) 








K.(ia,) = 


In any actual numerical computation the infinite product in Eq. (5.1) is, of course, to 
be replaced by a finite number of terms (recall that 8,,— mx as m — ~). The number 
of terms that is required to give an accurate answer is of course dependent upon / and 8, . 


TABLE 2. 








n Ba(h = 2) Br(h = 15) 

0 1.1828 +% .5832 1.5638 +7 .1046 
1 3.3106 +2 .6499 4.6886 +72 .3234 
2 6.3014 +7 .3277 7.8016 + 9749 
3 9.4248 +7 .2138 10.8705 +2 .9070 
4 13.7139 + 2 1.3629 
5 16.1332 + 2 1.4236 


18.9644 + 7 1.0572 
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In this paper g.(x) was computed for h = 2 and h = 15. Values of 8, for h = 2 and 


h 15 are given in Table 2. Let us first consider the case h = 2. Then the real parts 
of ay and a, are given by .59 and 1.65 respectively hence we may expect the fundamental 
term to give an accurate result for x > 2. Carrying out the necessary computations we 
obtain 

Pf n(x, t) = .61 exp (— .59z) exp [i@t — .885 — .292x)] ” 

; (5.2) 


+ 0 [exp (— a,2)]; h = 2. 
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Fic. 3. Ground-water phase lag vs. distance inland, z. (Dotted lines indicate the expected correction 


to the computed curves when all modes are considered) 


In Figs. 2 and 3 the amplitude and phase lag in the ground-water fluctuation have 
been plotted. The extrapolation of the results to = 0 are indicated by dotted lines. 
Actually it is not difficult to obtain another term in the series, but unless particular 
quantitative information is desired for this value of h it hardly seems necessary to do 
that. It might be mentioned that three terms were more than sufficient in evaluating 


the infinite product in Eq. (5.1). 


In the case that h = 15, the fundamental term can only be expected to be accurate 
for x > 7. We obtain 
Pod w(x, t) = 074 exp (— .1042x) exp [it — 1.68 — .007z)] (5.3) 


qi 


+ 0 fexp (— a@,2)]; h = 15. 
In order to obtain results valid for x = 1 or 2 when h = 15 would probably require the 
computation of three or four terms of the series. However in view of the results for 
h = 2 and these results for x > 7 it is clear that the amplitude curve for h = 15 will lie 
almost exactly on the curve given by the infinite depth theory” (see Fig. 2). 


*The amplitude and phase lag curves for h = © have been taken from the results given in [2]. 
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DIMENSIONLESS WAVE NUMBER 














Fic. 4. Dimensionless wave number vs. h; \ is the wave length predicted by shallow water theory; 
Xo is the wave length of the fundamental mode predicted by finite bottom theory. 


It is interesting to note that the shallow water theory and finite depth theory predict 
an exponential decay in x for the amplitude of the ground-water fluctuation; this is in 
contrast to the algebraic decay, (like x~'), predicted by the infinite depth theory. Also 
the phase lag predicted by the shallow water theory and finite depth theory continues 
to increase with x, while that predicted by the infinite depth theory approaches 7/2 
with increasing x. (This is clearly illustrated in Fig. 3). That g.(x) as given by Eq. (4.19) 
approaches the infinite depth result given in [2] as h — ~, cannot be seen easily from 
(4.19). However an examination of K(&) as given by Eq. (4.9) shows that as Ch — ~, 
K(t) 7 + (# + ze)'” which we might denote by K.(¢). This function is the one that 
occurs in [2]. It is interesting to note that K.(¢) has singularities of the branch point 
type, and in the limit as e — 0 these singularities will occur at the origin. This explains 
the algebraic behaviour of n(z, ¢) for h = ~. In contrast for any finite value of h the 
strip of analyticity of K(é) is finite even when e — 0, and its singularities are poles 
rather than branch points; hence the exponential sort of behaviour for n(z, ¢) for finite h. 

A plausible physical explanation for the fact that the amplitude curves for the 
ground-water fluctuation lie continuously below one another as h decreases (see Fig. 2) 
is the following. Imagine that our porous medium and fluid occupy the strip — H < 
y <0,-— © <2 < o. Suppose that we apply a uniform pressure on the half line y = 0, 
— © <2 <0; then fluid in the left half strip will be forced through the gap — h < y < 0 
and the free surface given originally by y = 0, x > 0 will rise. The amount of fluid that 
can be forced through this gap, and hence the effect that the pressure variation can have 
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on the free surface, is proportional to the gap distance, h. So with decreasing h the 
amplitude of the free surface fluctuation is lower and dies out more quickly. 
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BOOK REVIEWS 


Geometric algebra. By E. Artin. Interscience Publishers, Inc., New York and London, 
1957. x + 214 pp. $6.00. 


This incisive monograph leads one by rapid and brilliant steps through the author’s own algebraic 
foundation of affine and projective geometry, into a bird’s-eye view of orthogonal and “‘symplectic”’ 
geometries over general fields and, after a brief discussion of determinants over non-commutative 
fields, the structure of what modern algebraists now refer to as the “classical” groups associated with 
these geometries. The generality and penetration of the results obtained are truly remarkable. Thus, 


’ (generalized dot products) 


“fields” are not even assumed to be commutative; ‘‘metric structures’ 


need be neither symmetric (though x-y = 0 is assumed to imply y-x = 0) nor positive. Another inter- 


esting feature is the systematic use of coordinate-free methods. 

On the other hand, ‘“‘applied’’ mathematicians primarily interested in the real and complex number 
systems will probably find the book unnecessarily abstract and sophisticated, and the author’s antipathy 
to coordinates and matrices (p. 13) exaggerated. Therefore, it is perhaps well to recall that all mathe- 
matics is potentially 


original fo le force, among 


ipplicable, and that the volume being reviewed is at least to be admired as an 
the very best works of its kind. 
GARRETT BIRKHOF! 


On human communication. A review, a survey, and a criticism. By Colin Cherry. Pub- 
lished jointly by the Technology Press, M. I. T., Mass., John Wiley & Sons, Inc., 


New York, and Chapman and Hall, Ltd. London, 1957. xiv + 333 pp. $6.75. 


This is a very attractive and useful book. On pages 2 and 3, the author explains his purpose, in a 
statement which we regretfully abbreviate, but using his own words: “‘At the time of writing, the various 
aspects of communication, as they are studied under the different disciplines, by no means form a unified 
study; there is a common ground which shows promise of fertility, nothing more./ We shall attempt a 
review, a survey and a criticism of the study as it is being developed./ This book is introductory./ 
The book is written for that curious person, the ‘general reader’./ We are seeking to extract from/ 
semantics, psychology/ the common 


different sciences—linguistics, phonetics, communication theory, 
related concepts and ideas concerning communication.”’ 

The applied mathematician, being usually a man who likes to look over the fence, will derive both 
pleasure and profit from this very clearly written, conscientious, and unassuming survey. He will get 
a good idea of the work of a very interesting group of people, located at Harvard and M.I.T., experts 
in the various disciplines quoted above, who for the last fifteen yea 
effort, have built up this broad common field of communication. The present book is the first volume 
of a series, “Studies in Communication”, published by M.I.T. Press and Wiley, and the editors of the 
series, William N. Locke, Leo L. Beranek and Roman Jakobson, are prominent members of the group. 
Claude Shannon, who has done so much in the mathematical part of the field, has recently come to 
Cambridge. Prominent visitors from England have been A. Tustin, D. Gabor and the author of this 
book. Their work and that of more than a hundred people, from Adrian to Zipf, is discussed in the 
book and referenced in a marvelous bibliography of 367 titles. 

To the applied mathematician who, after reading On Human Communication, would like to go fur- 
ther into communication theory, we would recommend reading, in order of publication, the excellent 
books of Shannon and Weaver (Illinois Press, 1949) Woodward (Pergamon-McGraw-Hill, 1953) and 
Brillouin (Academic Press, 1956). Some knowledge of cryptography (a habit-forming pastime) is very 
helpful. 


rs, in a remarkable cooperative 
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ON THE PERIODIC SOLUTIONS OF THE FORCED 
OSCILLATOR EQUATION* 


BY 
R. M. ROSENBERG** 


University of Toledo 


Introduction. The phenomenon of “subharmonic vibrations” has been considered 
by many authors. The research on the subject goes back nearly one hundred years, 
beginning (probably) with Melde [14]}, Helmholtz [5], and Lord Rayleigh [17], all 
before the turn of the century, and continuing uninterruptedly to the present. For a 
fairly complete bibliography on this subject, the reader is referred to the classical paper 
on non-linear engineering problems by von Karman [8], and the books by Stoker [20], 
Minorski [15] and McLachlan [13]. ; 

In this paper, we shall discuss a single degree of freedom system whose mechanical 
model might be a mass under the action of an “elastic force” (linear or non-linear, 
restoring and/or exciting) and of a simple harmonic forcing function of frequency w. 
Moreover, we assume fhat small quantities of positive damping are present in the 
physical system, but absent from the equation of motion. Such damping can be shown, 
in general, to reduce the free vibrations to negligibly small amplitudes in a finite time. 
Consequently, the periodic solutions of this system are those usually referred to as 
“steady state vibrations.’’ We are assuming here that the effect of small damping on 
the motion is slight, but later on we shall prove this to be the case. 

It follows that the equation of motion of the system is the oscillator equation 


x’ + f(x) = Po cosat, (1) 


where P, and w are non-vanishing, finite, real constants, x is a dependent real variable 
and ¢ is the independent variable time. For the present, the function f(#) remains 
undefined. 

Equation (1) and its harmonic and subharmonic solutions have been explored widely 
and by a considerable variety of methods, but it seems that the treatment was restricted 
in most cases to almost linear, and odd, f(x). Ludeke [12] is one of the early investigators 
who has considered some cases, and experimented with some models, in which the 
departure of the elastic force from linearity is significant. However, he only considers 
very special cases; for instance, his treatment of the 1/3 order subharmonic vibration 
is restricted to special values of the forcing frequency. Moreover, his experimental model 
may be subject to an equation with periodic coefficients rather than one of the type (1) 
because the inertia of his pendulum varies harmonically with time. McLachlan [13] has 
given an example of a special equation of the type (1) in which f(z) is strongly non-linear, 
and this equation was also given by Cartwright and Littlewood [1]. Moreover, there 
is a great deal of significant literature [6, 7, 9, 11, 16] exploring solutions of equations 
like (1). In an earlier paper [19] the origin of subharmonics of odd orders was discussed; 
however, that investigation contained many conjectures, and the treatment was re- 
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stricted to odd orders. Here, we would like to lay a proper foundation for the earlier 
speculations, and to generalize the treatment to include the even-ordered subharmonics. 

We shall be concerned only with periodic solutions of (1), and we use the following 
definitions. (I) If the Fourier expansion of the periodic solution of (1) contains a term 
A, cos (w,t + ¢,), where w, = w/r and r is an integer, and if A, ¥ 0, the solution is said 
to be subharmonic of order 1/r. (II) If the solution is subharmonic of order 1/r it may 
be written as 

x = A. cos (2 t+ e.) a > s A; cos (w;t + ¢,). 
ine 
If in this solution |A, | >> | A, | for all 7, the solution is said to be a strong subharmonic 
of order 1/r. (III) If A, 0 and A; = 0 for all7, the solution is called a pure subharmonic. 
(IV) If A, # 0, A; = 0 for all 7, and gy, = 0 as well, the solution is said to be a simple 
subharmonic. 

In view of these definitions, the case of r = + 1 is not excluded. Therefore, the 
harmonic response is in the class of subharmonics; in fact, it is the subharmonic of order 
1. Moreover, it is quite possible for a solution to contain Fourier components of several 
distinct subharmonic frequencies. Such solutions will be called multiply subharmonic. 

Special interest centers on periodic solutions of (1) which arg of period 2xr/w where 
r is an integer. They are subharmonic in the sense of our definition and arise when 
A; = 0 for all i < r, except perhaps when 7 = 0. Such solutions may be written in the 
form 

z=B+ Dd B; cos (# 1 + ¢,] 
7=1,2,*** 
and they will be called subharmonic steady states. They are, in general, multiply sub- 
harmonic unless r is a prime number, and B, = 0. Evidently, the simple and pure 
subharmonics are subharmonic steady states. 


The simple subharmonics. Under the restriction that f(x) is analytic we will show 
that to every simple subharmonic 
= Le cos — ‘: (2) 


where 2, is a non-vanishing, real, finite but otherwise arbitrary constant, and r is an 
integer, there belongs one and only one equation 


x’ + f(x) = Py coswt (3) 

capable of producing it. Furthermore, the “elastic force’’ in that equation is of the form 
f@= Lv2" (4) 

and the coefficients y? , (n = 1, 2, --- , r) are uniquely determined by the integer r and 


by the parameters of the equation. Moreover, any “equation of a simple subharmonic” 
(i.e., any equation possessing a simple subharmonic as a solution) is reducible to a one- 
parameter equation. 

In view of (2) and a well-known trigonometric identity we have, when r > 2 is an 
integer*, 


*Excluding negative integers is no restriction since the cosine is an even function. 
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r r-2 _ 
cos wt = 2(Z) at 2-*(2) 41-3) (2) 
Lo 1! Le 2! Xo 


_%- 9-8 a(t) 4 thr — Sir — 6 — 7) y(t) - 


3! Lo 4! Xo 








with the series terminating when one of the coefficients vanishes. For r = 1 and r = 2 
we have, respectively, cos wt = 2/2» and cos wt = 2(2/z,)? — 1. Since 2° = — (w/r)*x 
we find for the elastic force, when r > 2, 


, row (x fz \’ r s(x \"? , rr — 3) (x \"* 
fe) = —- |—) + Pol 2) —- G2 ee | 
r \Zo Lo 1! Zo 2! Xo 


_ rr — 4)(r — 5) a(2)" 4 r(r — 5)(r — 6)(r — 7) a-(Z)" wee |, 


3! Xo 4! Zo 


(5) 








When r = 1, f:(z) = (ow + Po) (x/20) and when r = 2, 


filz) = za (=) + p, a(%) “ 1]. 


It will be seen that the coefficient of P, is in every case a polynomial of degree r con- 
taining only even, or only odd, powers according as r is an even or an odd integer. For 
every positive integer r, the elastic force is also a power series of degree r. When r is 
odd we have 


se) = 2 (Z)4p, © a2) 6) 


> —_— 
y Xo n=1,3,°¢* Xo 


and when r is even, 


2 r n 
Low (x al B wr 
y= (Sar S af®) , 
. r° Lo ~~ S Lo P ( 

where the a® depend on r only and are defined in an obvious manner through (5). It is 
an interesting fact that f,(x) is an odd power series when r is odd, but it is not even for 
even Tf. 

We have shown that the set of equations having simple subharmonic solutions of 
orders 1/r are 


— r ay 
xo t+ i +P, >> at(=) = P, cos wt, (r = 1,2, ---). (8) 

n 0 
If we set wl = r, x/z) = £ and P,/(xqw') = k, the simple subharmonics become — = 
cos r/r, (r = 1, 2, ---) and (8) becomes 


E +r%ttk Dak =keoss, (= 1,2,-->). (9) 


These equations contain only the single parameter k. In view of the definition of k and 
the restrictions on Py , 2 and w it is evident that all finite, non-zero values of k are 
admitted, but k = 0 is excluded. In Table 1 we tabulate the coefficients a* for the equa- 
tions of the simple subharmonics of orders 1/r, (r = 1, 2, --- , 9). 
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TABLE 1 


1 0 4 6 7 9 
a ar Q@ Qr a ar ar a ar 
] 0 l 0 0 0 0 0 0 0 0 
2 —] 0 2 0 0 0 0 0 0 0 
3 0 3 0 t 0 0 0 0 0 0 
$ 1 0 —8 0 8 0 0 0 0 0 
5 0 5 0 — 20 0 16 0 0 0 0 
¢ l 0 18 0 — 48 0 32 0 0 0 
7 0 7 0 56 0 —112 0 64 0) 0 
8 I 0 32 0 160 0 — 288 0 128 0 
9 0 9 0 -120 0 132 0 -—576 0 256 
The a-coefficients in Table 1 may be constructed from the following relations: 
a, 0 forall n>m 
0 when m+n = 2p — 1, (p = 1, 2, ---) 
a, = 0 when r= 2p — 1 
= —] when r = 4p — 2/(p = 1, 2, ---) 
= +] when r = 4p — 4) 
a, = r(—)""?”" when r = 2p — 1, (p = 1, 2, ---) 
a, = 2a’ ie , ip = 2, 3, ) 
an, = Qa, —a,.» when m>n+2 
when m+n # 2p — 1, (p = 1, 2, 

Finally, we observe that it is evident from the construction of the equation of the 
simple subharmonic £ = cos r/r, (where 7p is a positive integer) that, if f(€) is analytic, 
there exists one and only one equation of the form &* + f(~) = k cos r capable of pro- 
ducing it; this is the equation &* + f,,(&) = k cos r, where 


To 
f.. =r? +k > are" 


and the a”. are defined above. 
We would like to discuss briefly the significance of the equations of the simple sub- 


harmonics. 


If we are asked to find the “steady-state” periodic solutions (if any) of x° + 
>, o'2” = P, cos wt, we can do this by straightforward, elementary means only if 
; xl ax” = a'x. For that case we find x = x, coswt, where 17) = P,/(a’' — a’), [(a’)'” ¥ +e] 
Now, this last relation is an equation of the sort g(P, , a’, w, 2) = const., solved for 
the amplitude z, . The method of obtaining this important relation involves techniques 
which fail unless a? = a® = --- = O. In fact, the failure of these techniques in every 


but the linear case, presents one of the fundamental difficulties in the discussion of 
non-linear equations. 
° . _ o : : 1 
There is, however, another obvious way of finding the relation g(P> , @, #, %) = 
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const., and this way can be taken in the non-linear as well as the linear case. If we are 
asking for that function f,(x) for which x + (2) = P» cos wt has the solution x = 2 
cos wt we find by direct substitution that f,(z) = a'z, where a’ = (Py + w x) /Xo ; 
but this last is precisely the equation g = const. solved here for a’ instead of x, . If we 
make the transformation from ¢ to 7, from x to &, and put P,/(aw) = k, the desired 
solution becomes & COS T, if (&) = a’é, and a’ = j +. k. In view of the arbitrariness of 
ro , Po and w, this latter result preserves the generality of the former. As a consequence 
of this method of establishing a correspondence between equations of the form £° + f(£) = 


k cos r (with f(£) analytic) and solutions = cos r/r (with r and integer) we find that 
the linear, sinusoidally driven oscillator equation is not the only one having a “‘steady- 


state”’ sinusoidal response; instead, there is an infinite set of such equations. They are 
the “equations of the simple subharmonies,”’ and the linear equation is simply that of 
the simple subharmonic of order 1. 

It may be disappointing to those who cherish general solutions to equations with 
arbitrary parameters that the equations of the simple subharmonics have very special 
coefficients. We share this disappointment. However, a moment of reflection shows that 
these ‘‘special coefficients” actually define the relations which must exist between the 
parameters of the equations and the amplitude and frequency of the “steady-state” 
response. In this sense, these equations are no more special than the linear oscillator 
equation under simple harmonic excitation. Even the technique of prescribing amplitude 
and frequency of the solution before obtaining additional results is familiar in non-linear 
vibration problems. For instance, when searching for the harmonic response of the non- 
linear harmonically driven oscillator equation, the relation which exists between the 
amplitude of the response (of prescribed frequency!) on the one hand, and the frequency 
of the excitation on the other is normally obtained by prescribing the amplitude and 
computing the frequency [20]. Here, we have merely gone one step further and have 
inserted the governing relations in the equations of motion. 


Stability of the simple subharmonics. It is well known that the stability of the 
periodic solution £ = cos r/r of (8) depends on the stability of solutions of a Hill equation. 
In the case of the simple subharmonics, the stability investigation gains from discussing 
the odd and even subharmonics separately. Accordingly, we begin with an examination 
of the stability of the simple, odd subharmonics. 

The stability of the simple subharmonic £ = cos r/r of the equation 


r 


+rvt+k Dd) a" =keost, (r = 1,3, ---) 


n=1,3,°°° 


is identified with the stability of solutions of 
[_, : 
ui tr? tkir(—)oP? + D> na? cos” +/r | pu = 0. (10) 
\ n=3,5,**° 


In this equation, the term belonging to n = 1 has been written ahead of the summation, 
. -1)/2 
and use has been made of the relation a! = r(—)‘"~”” for all odd r. 
When m > 2 is an even integer, we may write 


m (m=-2)/2 
= T ~ im l ~'" (m < * 
cos" — = 2 of! m|J+ > ("") cos (m — 21) *} 


~ 9 #=0,1,2,¢*° 
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and when m > 3 is an odd integer we have 


(m—-1)/2 = 
es in ree m - 3 
cos" — = 2 + (™) cos (m — 2)-, 
Tr ° z ae 


t=0,1,2, 


where the conventional notation for the binomial coefficients has been employed. Applica- 
tion of the first of the above relations reduces the stability problem of the odd simple 


subharmonics to a discussion of the Hill equation 


| : . f | ieee = ane {1 n— ] 
wi + (4? + khr(—) + dh wil5(e-! 
' n=3.5,°°° &@ _ 2 
(n-3)/2 
re 2. bs ; ‘) cos (n — 1 — 22) |} = 0. (1) 
¢=0,1,2,°°° \ 


A quantitative discussion of the stability of this equation is not possible in view of the 
complexity of the periodic coefficient. However, there exists evidence [22] that under 
certain circumstances the stability of the solutions of (11) is not highly sensitive to 
the precise form of the periodic coefficient. Accordingly, we shall adopt the viewpoint 
that the coefficient of u is a periodic function which, for purposes of a stability analysis, 
can be adequately represented by the leading terms of its Fourier expansion. One of the 
circumstances which must be met is, obviously, that the leading term must have the 
largest coefficient. If we write the second sum in (11) in an ascending order of Fourier 


terms we have 
—" [n-1 a 
S ("—") sm —1 - 29 
POPS: oe GORD D r 


2 tr (n — 1)r 


l 
2T “ 
n—3]cos—+in-—5 cos — + +--+ + cos ——. 
sat r a 6 gna r rT 


In this sum, the magnitude of the coefficients decreases continually—i.e., 


in — 1 n—-1 
ia 2) eee 

. -* |\2z— ] 
for all integers iin 0 < 7 < (n — 8)/2. If only the leading term of the Fourier expansion 
of 


(n 


~/* n — 1\ ; Pe 
>. ( ) cos (n — 1 — 22) - 
1 j Tr 


i=0,1,2,°°° 
is retained, the stability of the odd simple subharmonics depends on that of the solutions 


of 


du : 


qd (aj + 65 coszju = O, 
where 
1 2 r ? n = 1 
br “a)/8 1a, 
a 4 ot 4 a aad ee - rm pent |? w - 
n=3 5 ** _ 
Y 
(12) 
, n{n—1 
. kr° na 
bi = ln — 3}. 
2 n=3,5 2 . — as 
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Now it can be shown that 


n—1 


‘ i _ F = _ 4) 
> a (r=) = r when r= 4m ‘ oe ee 


n=3,5,++* & 2 = 0 when r=4m-—3 


- nae (* ~ : 
= i n—3]=r. 
ee —— 


As a consequence, the coefficients in (12) become simply 
] ? ) 


tn FI ~ rb) 
bb = 51k 


The elimination of the parameter k between them yields for all odd integers r > 3 the 
relation 


5 5 ia 
% =Z+5 bo. (14) 


We call (14) the “stability characteristic” of the simple odd subharmonics of order 
1/r. It has the remarkable property of being independent of the order 1/r. This stability 
characteristic is shown together with the stable and unstable domains of the ab-plane 
in Fig. 1. From it, we see that simple odd subharmonics have a stable range when 
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1/4 < a3 < land 0 < bj < 3/2 where the upper bounds are somewhat larger than 
the least upper bounds. As aj and bj increase, there follows a fairly extensive unstable 
range for 1 < aj < 2 and 3/2 < bj < 5/2, but near aj = 2, bh = 5/2 a very narrow 
stable range exists. This is again followed by an extensive unstable range and an extremely 
narrow stable one near aj = 4, bj = 15/2. 

One of the interesting features of this result is that the simple odd subharmonics 
are never stable when k is negative. In fact the bounds on k which insure stability (in 
* where, again, the upper bound exceeds the 


the first stability range) are 0 < k < 3/r 
positive k implies a hard 


least upper bound somewhat. In the case r = 3, for instance, 
spring, negative k a soft spring;* it is seen then that simple subharmonics of order 1/3 
can not occur in a system with a soft spring. 

Another interesting result is that the simple subharmonic may be stable when the 
equilibrium position £ = 0 is not. For instance, in the case r = 3, the equilibrium position 
becomes unstable for k > 1/27 while the subharmonic is stable until k = 3/27, nearly. 

We examine, now the stability of the even simple subharmonics. When a procedure 
is followed exactly analogous to the preceeding one, except that r is even, it is found 
that the stability of the even simple subharmonics is identified with the stability of 


solution of the equation 


= + (a + 05 coszju = 0 


and the coefficients are 


as = I, Db, = 2r°k. (15) 


It is seen that the stability characteristic of the even simple subharmonics is a straight 
line parallel to the b-axis and crossing the a-axis at a = 1. From Fig. 1, where this stability 
characteristic is shown, we deduce that the even simple subharmonics are almost never 
stable except in some very narrow ranges the first of which occurs near b; 8 when 
k = + 4/r’, approximately. 

The strong subharmonics. In this section we inquire into the properties of solutions 
of equations which “lie in the neighborhood” of those of the simple subharmonics. 


Specifically, we consider the equation 


E g(7,é, [ule yu gir +T,é,|ul&,n), (16) 
with | uw | small and g(r, é, 0, 0) = —/f,(é) + k cos r. The function g is analytic in all 
variables. When yu 0, (16) has the solution = £&)(r) £,(r + T), but T need not be 
the least period of & or g. When yu does not vanish, (16) has the solution é n(r, uw). 

n(r, u) Which is periodic in 7 of 


It can, then, be shown [2] that (16) has a solution & 


period 7 and analytic in both arguments. Moreover, there is only one such solution for 


each u. It should be noted that (16) includes as a special case the equation £° + | u 
¢ + f,(&) = k cos 7, so that it has been shown that small positive damping does not 
have a profound effect on solutions obtained when damping is neglected. 


Since n(r, u) is analytic in the small parameter u one may write 


+ 


no(7, 0) + Do w'n(r, 0) = f(z) + Do w'é(7), 


In hard springs, the stiffness increases with deflection; in soft springs it decreases with deflection. 
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where every £; must be periodic since 7 is periodic. Thus, we have 


si 
: T , ir . jr = 
y(t, uw) = cos , + p Oe ps u'( a, cos 7” + 6,; sin ir) (17) 
t 7 
and, in view of the earlier definition, n(r, «) is in general a strong subharmonic steady 


state. 
As a special case of an odd strong subharmonic we consider a generalization of 


Duffing’s equation 


gs 4 pe 4 if ‘le art” + ut’*? = k cos 7, (7 = 1,3,---,), |] small. (18) 


It will be noticed that this reduces to Duffing’s equation when r = 1. 
When the theory of (17) is applied to (18), the solution to the latter, within linear 


terms in p, turns out to be 


a r + 2 r+2 . 
| ar+2 T . T 
S | —ywrt+1 /2 kr | cosa—+yp > a,; cos = , (19) 
“<2 r 7=3,5,°+° Tr 
where for j 3,5,°*:, <(n—1)/2 
r+2 \ : r n— 1 } 
; si . na, 
aq,=—-|r+2- -5)/ %3 Q—-jfyr?+hrt+ > gets — | : 
9 / / n=3,.5,0+* & 2 = J \ 
and for y] > (n — 1)/2,°--- ,r—2,7r,r+2 
r+3 
a; -(r+2-i)/ 2a - fr + kr]}, 
2 
where the usual notation for the binomial coefficients has been used. The solution (19) 
is in general — subharmonic. Suppose r = 1, to be odd integer which is not a 
prime. Let it have factors jo , J: , *** , Jm Such that the ratios ro/jo: , -** » To/ jm are 


intewers (yp , *** . Cy Then (19) becomes 


0 


r ' « / m ° 
, , T T : r 
g [ —plr +1 / 9 in| COs ~ + py _ Q,;; COS— + yp > 2 ay; cos = : (20) 
> , 0 jis 


‘<0 C; 


where the first term is strongly subharmonic of order 1/7, , the terms of the first sum are 
subharmonic of orders 1/e; , (¢ = 0, 1, --- , m) and the terms in the second sum are not 


subharmonic. This explains why one can find subharmonics of many orders in the 


solution of a non-linear oscillator equation. 
The stability of the solution (19) depends on that of a Hill equation. When only 


the leading terms of the periodic coefficient are retained, the “stability equation” 


becomes 


g vu 


s+ [(a5 + wal) + (65 + ub}) coszju = 0, 
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where a; and bj are given in (13), and 


r+2 7 
a; = E r+ i /2* |r +.1)[1 — (—)°°?7)/4 + & + 2)}, 








r+2 1 r+2 
bt = Ei r+1 |x 7 Kim — 4)/4+4r-—1]rt+ 2}, 
2 = 2 
r 2n n = 1 
~ Na, eer 
no = Zz. waar n— 1 / wt — (-)"°?”*}}, 
ase ene —— 
: no” n—1 
rn, = >. ni aia 3 /r. 
meS3 5le+e & 9 


It can be shown easily that n. and n, as well as aj and bj are positive for every r. 
The relation between aj and bj becomes a little more transparent when one writes 


r+2 r+2 
r+i}=|r—1t Jee) 
2 2 


in the expression for bj . Since 


r+2 r+ 2 
‘ei icir=—% 
2 \ 2 
the quantity «(r) > 1. In fact, «(3) = 2 and x(r) — 1 monotonically as r — ©. We now 
have 
r + 2 i 
b; r° r+ / 2 Wl — 4)/4 + x«(r)“(r + 2)} 
9 


from which it is seen by comparing it with the expression for aj that aj/bj; — 1 as r 
increases. For the likely conjecture that, nm» = n, = r/2 approximately, the relation 
between a’ and b] may be expected to be as shown in Fig. 2. It follows, that the stability 
characteristic of (18) lies everywhere in the neighborhood of the ‘‘odd stability character- 
istic’ of Fig. 1 but has a slightly different intercept with the a-axis. As a consequence, 
stable strong subharmonic steady states can occur for negative k. In particular, when 
r = 3, strong subharmonics may be stable in systems with soft springs while simple 
subharmonics are not, as shown earlier. 

When the case of even-ordered subharmonics, analogous to the generalized Duffing 


equation is examined the equation may be written as 


. 
ge trte+k(-) "+k Zz. art” + pt’** = k cos 7, (ry = 2, 4, +>) (21) 
guano 


and the solution within linear terms in » becomes 
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r’ r+2 r r+2 i. 
t= —poal(rt+2)+ecstt+u YL a, cos, (22) 
2 “— oO j=2,4,°° e 
where 
r+2 
a; = — r+ % -j / era — jyr*}}. 


The stability of this solution depends (within the earlier simplification) on the solutions of 


d* u 


de + [(@ + wai) + (65 + ub') cos zju = 0, 






































where @, and 6% are given in (15) and 
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The effect of these quantities on the stability can be established from the observation 
‘even stability 


‘ 


that d; ~ 0. Since | u | is small, the stability characteristic will lie near the 
characteristic’ of Fig. 1. However, it crosses the a-axis ata = 1 + bd, . Whether ua; 
is positive or negative, the characteristic will pass through a stable range of the ab-plane 
for some small! values of k in the interval > 0. Thus, while simple even subharmonics 
are almost never stable, strong even subharmonics are easily stable. 

Concluding remarks. We have concerned ourselves with the stable (and unstable) 
state of systems W host equations are special cases of 


~ F 


2 Me ot oa ig 


ctr 


é . F(é, & » by 


where F is analytic in all variables on any interval, but in general not linear in them. 
Trefftz has remarked [21] that stable solutions tend with t — © to periodic solutions 
of period r7 where r is an integer. In this paper, these solutions are called subharmonic 
for every 

In spite of significant progress toward a better understanding of this equation in its 
general form [6, 7, 9, 11], or when it is nearly linear [16, 18], one cannot in general predict 
which among a variety of possible periods the solution will exhibit. Nor is the 
‘‘mechanism” of subharmonic vibrations understood. With the mechanism we mean 
what Stoker [20] calls ‘‘a plausible physical explanation.’”’ An understanding of the 
mechanism seems to us essential for a profound and satisfying understanding of the 
phenomenon. 

The latter difficulty has been removed—at least in the case of the simple subharmonic 
steady states—by the formulation of the equations of the simple subharmonics. Evidently, 
if we are ready to believe that the equation of the simple subharmonic of order 1 is the 
only one capable of producing the simple harmonic response, we should be ready to 


accept a similar statement for the simple subharmonics of any order 


1 


Moreover, the formulation of these equations has robbed the linear equation of the 
singular position which it is usually thought to occupy as the only one whose solution is 
purely sinusoidal, or as the only one whose steady state has a uniquely determined 
period. For, it appears that the equations of the simple subharmonics of any order 1/r 
have a steady state & cos r/r of unique frequency 1/r. (We stipulate here that the 
steady state remains unchanged when the argument 7/r is replaced by r/r + 2zn, 
(n = 1,2, ---,r — 1);i.e., when the phase of the response is shifted by n cycles of the 
excitation.) This result is deduced as follows: the conventional method of determining 
the steady states (exclusive of their stability characteristics) is to assume as the solution 
a periodic function in 7 in the form of a Fourier series with undetermined coefficients. 
The coefficients are, then, fixed by the relations which are obtained when the solution 
is substituted in the equation [20]. If this technique is used on the equation of the simple 
subharmonic of order 1/r, and the period 7’ of the assumed solution is commensurate 
with, but larger than 2zr,—i.e., 7/2xr > 1 is a rational number—it will evidently turn 
out that all coefficients vanish except that of cos r/r whose coefficient is 1. When the 
period 7 is not commensurate with 2zr the solution is not periodic in 7 since the co- 
efficient of cos r/r does not vanish. It is, in fact, almost periodic. This result bears out 
the conjecture [3, 13], that subharmonics of order 1/r can occur in equations of the 


form &° + f(& i: cos r only if f(é) is a power series of degree r or higher. 


Finally, the multiplicity of steady states in non-linear equations of the type under 


discussion appears to have lost some of its mystery. An equation of the form £° + /(£) 
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k cos r either is, or it is not, that of a simple subharmonic of order 1/ro. If it is, the 
steady state appears to be unique. If it is not, the steady states are 1/r,,1/r,, «+: , where 
the r,; are integers and r; > r,_, . However, the steady states of orders 1/r, , 1/r2, ++: , 
are now merely the higher order terms of the Fourier development of a solution of period 


2rr, . This observation by itself is, however, not sufficient to explain the occurrence of 
a multiplicity of subharmonics. For instance, a solution>>, A, cos ir/r is periodic of 
order 1/m > 1/r only if m is a factor of r and, if the A; = O for all « = 
1, 2, +++ , (r/m) — 1. Now, it is obvious that distinct periodic solutions of one and the 


same equation have different initial conditions (since the equations considered here 

satisfy the conditions for existence and uniqueness of solutions when the initial con- 

ditions are prescribed.) Therefore, one can explain the existence of distinct periodic 

solutions by regarding the Fourier coefficients in the solution »¥ A, cos ir/r as functions 

of the initial conditions. Let the initial conditions of the 1/mth subharmonic be denoted 

by &(0), €(0) where m is a factor of r. Then the subharmonic of order 1/m will be 
I 


observed if it is stable, and if 


A,(é,(0), £.(0)) = +++ = Agymy—i(E,(0), &(0)) = 0 while A,,,, (&,(0), &.(0) ¥ 0. 


nm 


This observation does not aid in the construction of distinct subharmonic steady states; 
however, it provides a satisfactory mechanism for explaining them. 
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STABILITY LIMITS FOR A CLAMPED SPHERICAL 
SHELL SEGMENT UNDER UNIFORM PRESSURE* 


BY 
ROBERT R. ARCHER 
Massachusetts Institute of Technology 


Summary. An integration procedure for the differential equations for the finite 
deflections of clamped shallow spherical shells under uniform pressure is developed. 
Stability limits for the clamped shell are obtained for a range of the central height to 
thickness ratio from about 1 to 35. This serves to correct and extend previously known 
stability limits for this problem. 

1. Introduction. The existing literature on the subject of spherical shells allowing 
finite deflections under uniform radial pressure or point load at the apex divides into 
two parts. One part represented by the work of von Karman and Tsien [12, 13, 15];** 
Friedrichs [4]; Yoshimura and Uemura [14, 16]; Mushtari and Surkin [6]; and Feodosev 
[3] involves a determination of buckling pressures by means of a minimization of a potential 
energy expression for the shell with respect to a special class of deflection functions. 
Because of the rather special form of the assumed deflections in these papers, it is difficult 
to compare these results with integrations of the non-linear equations which as is noted 
in [8] can be derived as the Euler equations of the variational problem to minimize the 
potential energy of the shell; and therefore whose integrals correspond to a minimization 
with respect to a completely general class of deflection functions. 

The other part represented by the work of Biezeno [2], Kaplan and Fung [5], and 
Simons [9] is based on integrations of non-linear differential equations corresponding 
to those which are used in this paper for shallow spherical shells. However, since Biezeno 
integrated the equations after assuming special forms for the non-linear terms in the 
differential equations, it is difficult to decide what influence this has on the results. 

Kaplan and Fung are able to get integrals of the non-linear equations, but unfortu- 

nately they are able to determine buckling pressures, stresses, and deflections only for 
very low shells where the deflection shapes are of a simple type. In this range, their 
results are correct as far as they have gone in the perturbation of the non-linear equations, 
but appreciable corrections are to be found in the higher perturbations even in this 
range. 
Simons generalizes the power series method given by Way for flat plates (see [11], 
p. 338) to shallow spherical shells. Numerical results when compared with those of 
Kaplan and Fung for the clamped Shell are found to differ considerably owing to the 
retention of only a few terms in the power series solution. 

In this paper (and in [2, 5, 9]), the so-called “classical criterion” of buckling as 
distinguished from the “energy criterion” developed by von Karman and Tsien is 
applied to interpret the buckling phenomenon. In the “classical criterion”, it is assumed 
that a given state of equllibrium becomes unstable when there are equilibrium positions 
*Received October 26, 1956. The material in this paper was submitted as part of a thesis by the 
author to the Massachusetts Institute of Technology, May 1956, in partial fulfillment of the require- 
ments for the degree of Doctor of Philosophy. 

**Numbers in brackets refer to references at end of paper. 
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infinitesimally near to that state of equilibrium for the same external load. Thus, it is 
a question of obtaining the pressure-deflection relations for a given problem and properly 
interpreting the buckling pressure according to the above criterion. It is found that the 
center deflection to pressure relation used in [2, 5, 9] to interpret buckling must be 
generalized by interpreting buckling from a maximum deflection (in general, away from 
the center) to pressure relation in order to reveal the buckling in the cases where the 
deflection modes get more involved. 

It might be noted here that problems of finite axi-symmetric deflections of flat 
plates are included as a limiting case of the shallow shell, and thus the methods given 
in this paper carry over to these problems. 

With a view to the application of high speed digital computing equipment, the 
basic approach in this paper has been to reduce the integration of the non-linear differ- 
ential equations to the problem of solutions of algebraic equations by means of suitable 
sets of functions for the various cases. Thus, the rapidly increasing store of methods for 
applying computers to solving algebraic equations can be brought to bear on these 
problems. ‘ 

2. Equations for shallow spherical shells. The equations for the finite deflections 
of shallow spherical shells under uniform radial pressure which form the basis of this 
analysis are derived by Reissner in [8] and are listed here for reference: 


D/a(B"’ + B’/t — B/#) = —Y + Apa’t + B/E, (1a) 
1/Eha(v"’ + W’/é — W/€) = B — 46° /é, (1b) 
aV = ipaé, aq) —W + Spat + B/E, (2a) 
aN; = W/E, aN, = WV’, (2b) 
aM; = D(@’ + »B/é), aM, = D(B/é + vB’), (2c) 
u = E/Eh(W’ — v¥/é) w= —a [ Bdé, (2d) 
Eh° 

Vv = atH, D= (2e) 

si 1201 — 





Fic. 1. Element of the shell showing stress resultants and couples. 
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The equations of the middle surface of the spherical shell in its undeformed state 
are taken in the form 
To asin & Z = —acosé et. =e. 
where & is to represent one-half the opening angle of the shell. The components of 
surface loading for uniform radial pressure take the form 


p, = psing Pp, = —pcosd 
Equations (1) and (2) are the result of restricting attention to shallow shells where 
&. XK w/2. 


The following non-dimensional form of the variables will be used 


B* = B/é, , y* = W/Eht,/m’, (3a) 
p = —p/Per 5 x= &/,, (3b) 
N” = Em’a/h, m‘* = 12(1 — pv’), (3c) 
where p,, is the minimum buckling pressure for the corresponding complete sphere 
from the linear theory (p., = 4£h’/m’a’, see [10]). Using the non-dimensional variables 

(1) becomes 
1/r*L*B* + yp* = —2pr + y*6*/z, (4a) 
1/°L*y* — p* = —46**/z, (4b) 

where 

L*(--+) = (++-)" + (-++)’/a — (+++) /2”. (5) 


The corresponding expressions for stress resultants, stress couples, and displacements 


take the form 





aN, = Eh’ p*/xm’, aN, = Eh’y*'/m’, (6a) 
aM, = D(p*' + vB*/z), aM, = D(p*/x + v6*'), (6b) 
u/h = &./m7(p*’ — vy*/z)2, w/h = —dr m® [ 6 dz. (6c) 

z 

. ." 

| 

| 

: 











Fia. 2. Side view of element of shell in undeformed and deformed states. 
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3. Statement of the problem. The edge conditions for a clamped shell segment are 
given by 
u(é.) = 0, B(E.) = 0, (7) 


which in non-dimensional form becomes using (6c) 


y*’(1) — vy*(1) = 0, s*(1) = 0. (8) 
In addition to (8) for shells without central hole we have the condition 
B*(x), ¥*(x) regularat zx = 0. (9) 


The small finite deflections of a clamped spherical segment under uniform radial 
pressure are determined by solutions of the Eq. (4) subject to the conditions (8) and (9). 

4. Perturbation solution. A convenient method for getting solutions of (4) is to 
expand 6*, y*, and the inward pressure p into series in powers of a certain parameter 
and convert (4) into a sequence of systems of linear differential equations. The per- 
turbation parameter W will take the form of a ratio of deflection to thickness as a 
result of conditions imposed later. Thus, we write 


B* = > B,W', yt= > y,W', (10a) 
l=1 t=1 


ro) 


2p = > p,W' (10b) 


i=1 
and substitute into (4). Equating coefficients of powers of W to zero leads to the sequence 


of linear systems 


1/°L*8, + y, = —piz, (11.1a) 
1/N*L*y, — 6, = 0, (11.1b) 
1/\"L*B, + y2 = —p.x + Biv,/z, (11.2a) 
1/L* yy. — Be = —436i/z, (11.2b) 
1/L*8, + v= —pixt >> B.y;,/z, ene See ee (11-1a) 
1/\*L*y, — B, = —3 >> B,B,/z. (11-lb) 
i+j=] 
The boundary conditions in terms of 8, and y, become 
B,(x), ¥.i(x) regularat x = 0 (12.la) 
B81) = 0, Wi(l) — vy,(1) = 0 (12.lb) 


for all l. 
4.1. Determination of 6, ,¥,, and p,. 8, and y, are to be solutions of (11.1) satisfy- 


ing the conditions (12.1). We seek solutions of (11.1) as expansions in terms of Bessel 
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functions of the first kind in the form 


Bi(z) = Do an? Ti(rn2), (13a) 
n=1 
v(x) = C2 + pe by? Ts (And), (13b) 
n=1 
where the i, are defined by 
JX.) =0 =Wm=_1,2,---). (14) 


Thus, 8, satisfies (12.1), and it remains to choose c, so that y, satisfies (12.1). From 
(13b), we have 


Vi) =e + DL dToA) 


n=1 

and 

vyi(1) = ve, , 
where we have used (14) and the formula 

Fiat) = daoOhat) — 1/2d,042). (15) 

Now y;, will satisfy (12.1) if 

Ci = ao pe ade To(An)« (16) 

roa n=1 


If (13) is substituted into (11.1), it follows that 


—(A,/d)’an” +b,” = —C, + pdr. , (17a) 
a,” + (A,/d)*b,” = 0, (17b) 
where we have used 
L*[J,(Ant)] = —AnJi(Anr) (18) 
and equated coefficients of J,(A,x) to zero. The I, are defined by 
2 =D Pwdilut) (19) 
and using 
i ‘PI de = -1NSOd (20a) 
[ 20042) de = 4730) (20b) 
we get 
r, = -3 (21) 


~ And o(An) 
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Solving (17) gives 


' —2(c, + PA, d)° 


= - . —— (22a) 
XJ o(A.) [1 + (A,/d)*] 


We. ) ee 
2, + Pi) _ (22b) 


bo” = 
’ And o(An){L + (An/A)*] 


The coefficients p, will be determined from the following conditions on the deflection. 
Using the expression for w/h from (6c), it follows that 


w(x*) — w(1) . fr ee =e 
=-";/[ (6W+6,W? + ---) de, (23) 
h m Ji 
where x* (0 < x* < 1) will be taken at the point where the deflection is a maximum. 
If we impose the conditions 
ol 
1 = \°*/m | B,(a) dx, (24.1) 
QO = | B,(x) dz, (24.2) 
= | B,(x) dx, (24-1) 
then (23) reduce to 
‘ 7% —_ sf 
w(a*) w(1) W (23) 
h 


which serves to define the parameter W. The Eq. (24) will determine the coefficients 
p, in the expansion for the inward pressure p. 
To determine p; , we use (13) to write (24.1) in the form 


1 
. 


l ?/m? | > Gg JSilA Ft) Oe 


n=1 


_ 2, + pir > (An/A)*[Fo(Anw ) — Fo(An)] 


m os NT (An) [1 + (A,/d)*] 
or 
( = a ) 
— 2c, + D,) ms » rae eH , (26) 
where we have used the integral formula 
/ J(A.2) dx = - | J o(Ant). (27) 
(An . 


Since c, is known from (16), p, is determined. 
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4.2. Solution for any 6, , y, , and p, . Let 


B(x) = >» an? I i(Ant), (28a) 
n=1 
Vix) = ea + >. bo? J,0\,2), (28b) 
n=1 


then substitution into (11./) leads to 
— (A,/d)*'a,” + b,° = —(, + pT, + AL”, (29a) 
a, + (A,/d)*b, = GL”, (29b) 
where the G{’’, H‘” are defined by 


Gn =4D Ler DY ai?a;”), (30a) 


wicks (s,t = 1,2,--- ,l—1) 


BO an F ont 2 ONS. re (30b) 
~ l s+t=l 


#=1 j=1 


Il 


and the C:’ defined by 
1/rJ,(A,a)J,0\;2) = >> CiJ,0,2). (31) 


The G{” and HH” are known from the 1, 2, --- , / — 1 stages of the computation. Solving 


(29), gives 


I ( 2c, a Di)(An ‘)? (“) (1) y(L) /\\4]-1 
a,, = ) — DJ o(d.) —_ d H, 4. G, (1 + (A,,/A) ] (32a) 


pry _ 2 2. + pr) 


: 2a 4)-1 
LdaJo(a,) 1 Ha + On/d)'Gs bi + (r,/d)*] (32b) 


and for the /th stage (16) becomes 


l saat a 
Cc; al ¥ : } is Anda J (An) (33) 
—Va 1 
which forces y, to satisfy the boundary condition (12./). The condition (24.1) leads to 
an equation for the /th stage similar to (26) for the first in the form 
( 


| & Jaw =~ 30> 1" 
— We. »24 pO" a \; 
Aer + pi) =) De To(An)[1 + On a 





LS AAA? — GE J oQ.2") — Lime 
{> Mall + 04/04 : = 


Therefore, c, and p, are determined; and the a{" and b‘” are completely determined. 


Starting with aS”, b{”, ¢, , and p, ; the G\”’ and H°” can be computed and then 


) 


a, be, , and p. . Similarly, GS" and H‘” are computed leading to a{"’, b{”, c, , and 
p, . This can be continued to obtain any number of terms in the series for 8*, ¥*, and 
2p. Stresses and displacements can be computed from the Eq. (6). The buckling pressure 
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is determined from 
2p = > p,W' (35) 
t= 
by the condition 


dp/dW = 0. (36) 


5. Numerical solutions. This section involves numerical calculations using the 
integration procedure set up in the previous section to obtain new results for the finite 
deflections of a clamped spherical shell segment under uniform pressure. The perturb- 
ation of the non-linear equations is carried out far enough to determine the deflection 
curves as a function of the pressure for maximum deflections up to about one thickness 
of the shell. For the case of inward pressure, the buckling phenomenon is observed in 
this range of deflections and buckling pressures are found using the “‘classical criterion”’ 
for a range of the central height to thickness ratio from about 1 to 35. (See Figs. 3 and 4). 

In the numerical computation for this problem or for other cases of loading and 
edge restraint involving shallow spherical shells, it is necessary to compute the expansion 
coefficients, c,’, which enter the computation from the non-linear terms. These were 


computed directly from the definition 














1 
hin = BvIOT" | J (Acar) J (Ajax) J (An) dx (37) 
0 
1.2 
( So 
0.8F 
p 6 ° 
04+ .- - 
Oo 4 n 4 4 4 4 4 
2 4 6 = 10 iS 20 30 
d 
4 1 1 4 4 i L 1 i rl 4 4 
2 4 6 8 10 = 5 20 30 40 5 100 
H 
he 


Fic. 3. Known theoretical and experimental results are compared. A represents the theory in [5], B 
the theory in this paper, C the theory in [12], the solid dots the experiments in [5], the hollow dots 
the experiments in [12], and D the classical buckling load for a complete sphere given here for reference. 


using Simpson’s rule with one hundred values of the integrand. The results obtained 
by machine computation are given in Table 1. It is shown in [1] that these coefficients 
can be used to solve a variety of other problems involving other conditions of loading 
and edge restraint. 

In determining the buckling pressures given in Fig. 3, terms were computed in the 
perturbation series used in finding the pressure-maximum deflection curves until it was 
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seen that the higher order terms were no longer substantially affecting the curves in the 


neighborhood of the point where dp/dW = 0. The linear solution checks the linear 
solution given in [5] which is the same as the one given by Reissner in [7]. The second 


0.8F 





0.2- 





'@) i 1 1 1 rn 
| 2 3 3 5 6 
d 
Fic. 4. The curves A, B, C, and D represent results using two, three, four, and five terms respectively 
in the perturbation series of this paper. E and the cross are computed in [5] using two and four terms 
respectively in a perturbation series while the solid dots represent the experimental results reported 
in [5]. F represents the buckling pressure for a complete sphere (linear theory) and is given here for 





reference. 


TABLE la. Numerical values for the expansion coefficients C;’ and C3’. 


Ci’ X10 1 2 3 4 5 6 7 8 
1 10.123 3.0232 —0.20513 0.067211 —0.029812 0.015362 —0.0094689 0.0060496 
2 5.9445 1.9771 —0.15236 0.054438 —0.025949 0.014360 —0.0089100 
3 4.1592 1.4594 —0.11913 0.044558 —0.029812 0.012264 
4 3.1925 1.1564 —0.097467 0.037596 —0.032196 
5 2.5891 0.95663 —0.082360 0.032196 
6 2.1770 0.81565 —0.071173 
7 1.8778 0.71087 
s 1.6506 
C,’ X 10 ‘ 
l 5.4449 10.706 3.5608 —0.27440 0.098043 —0.046734 0.025862 —0.016047 
2 8.1792 7.4106 2.4816 —0.21304 0.082514 —0.042530 0.023878 
3 6.1991 5.5398 1.8982 —0.17083 0.069546 —0.036298 
4 4.88463 4.4079 1.5366 —0.14285 0.058814 
‘5 4.0093 3.6553 1.2906 —0.12160 
6 3.3932 3.1208 1.1123 
7 2.9386 2.7218 
2.5901 
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TABLE 1b. Numerical values for the expansion coefficients C3’ and C,’. 


(Vol. 


ci}x10 1 2 3 4 i5 6 7 
1 —0.53368 5.1437 10.821 3.7968 —0.30992 0.11592 —0.05795 
2 10.705 8.9547 8.0023 2.7420 —(). 24677 0.10046 
3 9.0110 7.3074 6.2026 2.1426 —0.20310 
4 7.2969 5.9928 5.0455 1.7617 
5 6.0597 5.0423 4.2489 

6 5.1395 4.3403 

7 4.4845 

8 

ci! x 10 

1 0.22866 —0.51834 4.9651 10.861 3.9340 —0.33159 0.12791 
2 4.6872 10.465 9.2271 8.3266 2.9027 —0. 26985 
3 9.5559 9.5422 7.8368 6.5981 2.3038 

4 8.4625 8.0664 6.5892 5.4518 

5 7.2569 6.7884 5.6391 

6 6.2774 5.9621 

7 5.5059 

8 

Taste 1c. Numerical values for the expansion coefficients C;’ and Ce’. 

Cs’ X 10+! 1 2 3 4 i5 6 7 

] —0. 12529 0.22878 —0.50064 4.8597 10.881 4.0204 —0.34613 
2 —0.49712 4.4294 10.286 9.3556 8.5295 3.0115 

3 10.020 9.6808 9.7889 8.1454 6.8637 

4 9.9644 8.9645 8.4968 6.9660 

15 8.8285 7.9038 7.3757 

6 7.7507 6.9680 

7 6.8521 

8 

Cr x 10" 

1 0.07686 —0.12983 0.2293 —0.48764 4.7862 10.892 4.0808 

2 0.22922 —0.47456 4. 2687 10.154 9 4264 8.6695 

3 4.1205 9.7031 9.6969 9 9222 8.3468 

4 9.6901 10.115 9.2316 8.7678 

5 9.4093 9.2271 8.2952 

6 8.5468 8.2610 

7 7.6093 

8 

TABLE 1d. Numerical values of the expansion coefficients Cy’ and C3’. 

Cc’ x 10! «1 2 3 4 i5 6 7 

1 —0.054954 0.083337 —0.12927 0.21819 —0.47798 4.7337 10.898 
2 —0.13705 0.22411 —0.46033 4.1587 10.057 9.4692 
3 —0.45308 3.9300 9.4784 9.6822 10.004 
$ 9.3001 9.6196 10.171 9.3924 
5 10.185 9.6224 9 4624 
6 9.5827 8.9196 
7 8.7609 





No. 4 


AN, 


~ 
0.031905 

—(0).052433 
0.085728 

—Q. 16981 
1.4948 
3.6684 
3.8044 
3.9618 


-0.063310 
0.11110 
22206 
1.9135 
1.6415 
4.9127 
5.2500 
1.8913 


8 
0.13531 
—0. 28374 
2.4148 
5.7337 
6.0342 
6.4774 
6.1939 
6.1189 


0.35609 
3.0890 
7.0547 
7.2246 
7.412 

7.4116 
7.4078 
6.9348 


Ss 
4.1255 
8.7706 
8.4868 
8.9559 
8.5534 
8.5929 
8.1406 


7.9766 
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C;’ 10+! 

] 0.039952 —0.058843 0.080990 —0.12290 0.21263 —0.47003 4.6946 10.901 
2 0.087558 —0.13310 0.21567 —0.44588 4.07865 9.9804 9.4976 
3 0.21762 —0.43106 3.7946 9.3122 9.6574 10.057 
4 3.7145 9.0100 9.5364 10.191 9.4949 
5 9.4823 10.178 9.7333 9.6154 
6 9.7833 9.7782 9.1540 
7 9.2635 9.0768 
8 8.5853 


order solution checks with that given in [5] for 3 < \ < 5 where the maximum deflection 
is at the center. 

The Massachusetts Institute of Technology digital computer, Whirlwind I, was used 
extensively to reduce the computation time. The programs used for the machine were 
checked by independent calculations using a desk calculator. 

In Figs. 3 and 4, the stability limits found in this paper are compared with known 
experimental and theoretical results. It is seen that the theoretical curve given in [5] 
based on two terms of a perturbation series is subject to considerable correction when 
more terms in the series are computed. In particular, the minimum value of H*/h for 
which buckling occurs must be revised upward from the value of about 0.67 given in 
[5] to 2.2 found in this paper. Also, it is seen that the experimental results depart from 
the theoretical results of this investigation as H*/h increases. As indicated in the refer- 
ences, this is probably due to the fact that the shell, under disturbances during the 
testing, jumps to a nearby buckles state before reaching the buckling pressure predicted 
by the theory in this paper which allows only continuous load-deflection processes. 
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ON TRANSVERSE VIBRATIONS OF SHALLOW SPHERICAL SHELLS* 


BY 
MILLARD W. JOHNSON AND ERIC REISSNER 
Massachusetts Institute of Technology 


1. Introduction. The present paper is concerned with the frequencies of free vibra- 
tions of shallow spherical shells of constant thickness. It has been shown earlier that 
for vibrations of shallow shells which are primarily transverse a considerable simpli- 
fication of the problem can be effected by a justified neglect of longitudinal inertia in 
comparison with transverse inertia [1]. 

Previous applications of this observation included a study of axi-symmetrical vibra- 
tions of spherical shells [2], and a study of inextensional vibrations of general shallow 
shells [3]. In the present paper we investigate vibrations of shallow spherical shells, 
without axial symmetry. Appropriate solutions of the differential equations are obtained 
and these are used to obtain the frequencies of free vibrations of a spherical shell segment 
(or cap) with free edges, in their dependence on the curvature of the segment and on the 
number of nodal circles and diameters. 

We may summarize certain qualitative aspects of our results as follows. Let H be 
the height of the apex of the spherical cap above the edge plane of the cap and let h be 
the wall thickness of the shell. When H/h = 0 we have Kirchhoff’s results for the flat 
plate. When H/h tends to infinity the frequencies of free vibrations of the cap tend 
either to a limiting frequency which may be called membrane frequency or they tend 
to the frequencies of inextensional vibrations which we have previously considered [3]. 
The membrane frequency is limiting frequency for all vibrations with one or more 
nodal circles. Its value is independent of the number of nodal radii and of nodal circles 
provided the latter is not zero. On the other hand, the frequencies of inextensional 
vibrations are a function of the number of nodal radii and presuppose that the number 
of nodal circles is zero. 

2. Differential equations for free transverse vibrations of shallow spherical shells. 
The differential equations which are to be solved are of the form 


VV'F — £ V*w = 0, (2.1) 
2 
DV?V7w + z VF + ph ae ~é. (2.2) 


The various quantities occurring in (2.1) and (2.2) have the following significance 
w = transverse (axial) displacement, 

Airy’s stress function, 

radius of middle surface of shell, 

h = wall thickness of shell, 

Eh, longitudinal stiffness factor, 

modulus of elasticity, 


by 
uu 


Si 
nt 
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368 
D = Eh'*/12 (1 — v’), bending stiffness factor, 
y = Poisson’s ratio, 
p = density of shell material, 
V? =()etr” ). + 7r°* ( )es , Laplace operator in polar coordinates r, 6. 
Stress resultants and couples, in polar coordinate form, are given by the following 
expressions. 
——.. gta if= 
ae or’ . r or r oo’ | = 
y (2.3) 
N - 0 (? oF) 
re ar \r 00 
- oV ~w 1 aM, Pres 
V,= —pD~ "4-94 _Ty,| 
or r 00 R . 
by (2.4) 
. oV -w oM,, “ae 
V, = —-D—— + —"- SQN, | 
roo or ki J 
‘O’w vy ow vy Ow 
M, _ - (2 a ee 3 «+ 4 
or r or rors | 
Vl D(? Ow 4 1 o’w + e). 9 &) 
My = tee Be ers (2.9 
, r or r 06 or | ) 
d (1 dw’ 
M,, = —(1 —»)D—(-) 
or \r o@ 


Equations (2.1) to (2.5), except for the inertia term phw,, , are the same as those for 
| 
\ 


4]. 


problems of statics 


3. Solutions for simple harmonic motion. We consider solutions of the form 


. cos n@ of ip ‘ 
(w, F) = [w,(r), f(r) ]e"”’. (3.1) 
sin n6 
Substitution of (3.1) into the differential equations (2.1) and (2.2) leads to the ordinary 


differential equations 


0 

LiF, — R L,w, = 0, (3.2) 

‘ . 1 , a 

DL,w, — php’w, + = LF, = 0, (3.3) 

R 
where the operator L, is defined by 
d° 5 d n- ‘ 
L, - dr” ¥ r dr 7 r (3.4) 
We may write the solutions of (3.2) and (3.3) as follows 

l ie 

Ww, = IMe T XM, (3.5) 


7 RDx 
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C C 
F, = ( + non) vt + Rat LiXe + be (3.6) 


In Eqs. (3.5) and (3.6) we have 


1 / E C0 1/4 a 
A = E (hop* _ {)| (3.7) 


und 
Ges C... logr, n= 09 
jen oe See (3.8) 
C7 oe, l<n 
C's.0 + Ca.o logr n=0 
‘i 3,0 T Uso 1027, : (3.9) 


Coa +t > l<n 


12C,.0r° + 4C. or (logr — 1), n=0 


lay 3 ' 
} aC; , + 4 oir log r, n=1 





yi = (3.10) 
 aikdi alg 
ee ee — a< 
"Find + “™4yq—n) =* 
x Ay nd (Ar) + Aon Y,(Ar) + Ag.ntn(Ar) + Ag. nK, (ar). (3.11) 


The functions J, , Y, , J, and K,, are Bessel functions and modified Bessel functions of 
order n. We note that the argument of these functions is real, as long as the frequency 
p is greater than a reference frequency which we denote by p.. and which is given by 


C 


i. =. (3.12) 
Po = hohe : 
Alternatively, p.. may be written in terms of the rise H of the shell segment which, from 
Z */2R, follows in the form H = a’/2R, where ais the base radius of the shell segment, 
as 
E\'? H , 
Do = a= 5° (3.12’) 
p a 


The solutions (3.5) and (3.6) involve eight arbitrary constants for each value of n. 
However, for n = 0 and n = 1 only six of these have physical significance. The quantities 
C;.5 and C;,, may be omitted as they do not enter into expressions for displacements 
and stresses. Furthermore, in order that the displacements u and v in radial and cir- 
cumferential direction be single valued we must have (see [4] for the case n = 0) the 


relations 
Cs.0 = 0, C3., = 0. (3.13) 
We next write resultants and couples in a form which corresponds to (3.1), viz., 


g ) 
1 COS N86) ine 
< e 


N 
ea 


= * (3.14) 
lsin n@} 
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In this way we shall have in terms of the functions F,, and w,, 


, : ue 2) 3 Fr. 
7 a7 7 (3.15) 
‘ dF, _ain 5 
bn = aT; Nw. = -—F, 
Now dr- . ¥ a . f 
; adL,,w n? d (w, r 
V.. _p| Ses _ , -; (~ N. 
i ! dr > ar\r | R Ven eee 
\ (3.16) 
. n d? (w,) ri. | 
— D| == w, F (1 — ») rat — = Nyon! 
D| =" u vn a3 \— | R Vre.n| 
_——s ew rm y dw, a mn " 
dr r ar r a 
dw d-w 
~ ee E —* — w+ ve | f (3.17) 
r ar 7 ay 
1 (w 
} = - — | — 
M.. +(1 )Dn © (2) 


4. Boundary conditions for spherical segment with free edge. In order that the edge 
r = a be a free edge the following four conditions must be satisfied 


N.. ta) = 6, Weta) = 0, (4.1) 
M,.,(a) = 0, y .. te) xe: (4.2) 


In terms of F,, and w, and indicating differentiation with respect to r by primes, these 


conditions assume the following form 


aF‘(a) — n F(a) = 0 : 
’ (4.3) ° 
nlaF/(a) — F,(a)| = 0 
a-w?'(a) +- qvwi(a) — vn w a) = 0 
. (4.4) 
a’ [L,,w,(a) |’ —ni(l — v) [aw;(a) —w (a)] - () 
Equations (4.3) may be simplified further, to read 
F(a) = 0, aFi(a) = F(a ; 
be (4.5) 
Fi(a) = 0, F(a) = 0: Zin 
In addition to boundary conditions for r = a we have regularity conditions for r = 0. 





Inspection of (3.5) to (3.11) reveals that these regularity conditions require the vanishing 
of four of the eight constants of integration in (3.5) and (3.6), namely 
Co, =(C,,= Arn = An = 0. (4.6) 


This leaves us with the following expressions for w, and F, 


l " y / oy 
= RD C. m a A, ad a(AP) + A 7 nt,(Ar), (4.7) 
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, n+2 Y 
F, (1 oe piena)C . a + C;..7° + me [— Ay nd (Ar) + Az nl n(Ar)]. (4.8) 

An equation determining frequencies of free vibrations is now obtained by sub- 
stituting the solutions (4.7) and (4.8) in the boundary conditions (4.4) and (4.5). In 
doing this it is convenient to consider the cases n = 0, 1 separately from the general 
and more complicated case 2 < n. 

5. Rotationally symmetric vibrations and vibrations with one nodal diameter. 
Substitution of (4.7) and (4.8) into (4.4) and (4.5) for n = 0 and n = 1 leads to the 
following system of simultaneous homogeneous equations for the constants of integration 
in (4.7) and (4.8) 


A, oludo(u) + 1 - v) J 5(u) | — A, oul o(u) —-(1— v) I 5(u)] = ( 


A, ,oF8(u) — Asoli(u) = 0! (5.1) 
| 


Ca 4 , ’ _ 
[ + at | 1,0 +- R [— Ay od o(u) + As ol b(u)] = 0 


9 eu 


) 


{ (j — v)( J (u) — pd {(u)) — p J(u) ] 
+ As .[(1 — »)(Ui(u) — wl i(u)) + w'l,(u)] = 0 





A, [C1 — »)(—J,(u) + a {(u)) + uJ 1(u)] 
' 5.2 
+ A; ,[(l — »)(—1,(u) + wl i(u)) — wu] = 0 (5.2) 
a 4 y 
l Ca ( 
& . ( 1 ~ 4 ) _— a , \ 
a ae | iT aRn [Ay (J i(u ped {(u)) | 
+ As alla) + ul) = 0) 
In these equations the quantity uw is defined as 
hpa* , Cor" ‘ 
=a =/- 7 — —; 3 
ila | D ? ~ DR? we 
or, with the reference frequency p.. of Eq. (3.12) 
hpa' 9 2) di (5.3) 
ro — : 5.3 
. p Pe 


It is noted that the systems (5.1) and (5.2) are such that in each case in order to 
determine admissible values of u it is sufficient to consider the first two equations of 
the system. These first two equations are the same as the corresponding equations for 
vibrations of flat plates and accordingly the equation determining possible values of yu 
is the same as Kirchhoff’s frequency equation for vibrations of circular plates with 


free edge, with zero or one nodal diameter, 


pu J(u) T,(u) | = 
~ ie | a fey n=0,1. 5.4 
2 | ee Th " ii 
We note from (5.3’) that values of p larger and smaller than p.. are accounted for by 


considering the following ranges of u. 


P>Ppo:p= ge"; r> ¢ m= 0,1,2,--- (5.5) 


’ 
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Dp = Dot b= — ns : > @, m= 1,3,'5,.*** (5.6) 


However, the characteristic equation (5.4) can be shown to possess no solutions of the 
form (5.6). Also, the substitution of (5.5) into (5.4) leaves (5.4) unchanged so that it 
suffices to consider only real positive values of u. We conclude that for n = 0, 1, the fre- 
quency p.. represents the lowest possible frequency of free vibration. 

For given v Eq. (5.4) is satisfied by an infinite sequence of positive values un = yu, , 
k = 1, 2,3, --- , where k is the number of nodal circles. Let p,,, be the corresponding 
frequencies of free vibrations and p<’, the values of these frequencies for the case of flat 


plates, for which p 0. According to (5.3’) we have 
- » { D ) , oye 
Paik = Mn [24 (5.7) 
and 
Pak = Pat + Do (5.8) 


We may write (5.8) in the form 


Pot = E + (25,) | (5.9) 
Dn,k Pn.k 


or with 
ite) Ee 2. Mey 
Dn k 7 pa D uy age 2 a 
(5.10) 
my = [1+ 80-4)" 
Dn.k L Ln k 
Known values of u, , , for vy = 1/3 and for k = 1, 2, 3, are given in Table 1. Values 


of the frequency ratio p,,.,/p.., as a function of H/h may be found in Figs. 2 and 3. 
These results are new only for n = 1, the case n = 0 having previously been considered 
in [2]. We note that the modes with no nodal circles (k = 0) correspond to rigid-body 
translations and rotations. 

We finally note that the frequency p,,, is equal to p.. in the limiting case D = 0, or 
H/h = ~. In this sense, the frequency p. may be designated as the membrane vibration 
frequency of the shell. 

6. Vibrations with two or more nodal radii. Introduction of the solutions (4.7) and 
(4.8) into the boundary conditions (4.4) and (4.5) leads, when 2 < n, to the following 


system of simultaneous equations. 


(1 +) ——— 8, + Ch. + a [-Ard) + Asal] =0, 6.1) 
\ uw / 4(n + 1) Lt 
{ x \ n+2 7 l a 
pS) OF, + Cha += [—Ar diy) + Asal] = 0, (6.2 
(1 tea in tm ia tol tnd (a) + As .ntn(u)] (6.2) 
- 
(1 — y)n(n — 1) CR. + {Il — en? — w?JJn(u) — (1 — udu) Aan ™ 
Be (6.3) 


+ {[((1 — vn? + pw’ l(u) — (1 — oul i(u)} As,. = 0, 
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Fic. 1. Curves » = p(x*) and » = 7 (x) when vy = 1/3, according to Eqs. (6.7) and (6.13). 

















H/h 


Fic. 2. Values of ratio of shell frequency pp,, to plate frequency p{°} in dependence on ratio of shell 


rise H to shell thickness h for various values of the number 2n of nodal radii and the number k of 
nodal circle s, for a value vy = 1 /3 of Poisson’s ratio. 
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Pak 
p (0) 
nk 
1,08 
1.04 
I. 
H/h 
Fia. 3. Enlarged version of part of Fig. 2. 
4 
ae K , ; ; en ; 
(l—vyn(l—n 7CUr. tT (11 — pn? J,(u) — [CL — vn + J(u) As (6.4 
».4) 


+ {1 — pn lu) — (1 — yn — uw |Ui(u)} As. 0. 


The quantities C* ,, and « are defined as follows. 


‘; ce ai cae 
1 oy a (6 


C = = deme Ce SS ars Cz, 5) 
Ra" Ra , 

and 
_ Ca’ , 18(] 2 [= (6 6 


The vanishing of the determinant of the system (6.1) to (6.4) leads to the following 


frequency equation 


‘g S,(u) ” 
p <p:5 mr — |, (6.7) 
k R,(u 
where 
S,(u) = 4n n? — 1) — vuln) liu) — Ji (u)T,(u) | 
(6.3) 
Ss n a n : 
+(n+1)01 - | 120) —_- 1.0) || 120) -_-- 1.0) | 
u yu 
and 
R,(u) = {1 — [uJ fu) — nm? J(u] + wd (uf td — vn" [ul )(u) — [,(u)) — w liu} 
— {1- yn" [ud 2 (u) — J,(u)] + pd i(u)} f(1 — v)[uli(u) — n’I,(u)] — uw T,(u)}- (6.9) 


' 
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In order to account for values of p both larger and smaller than p.. we consider the 
ranges of » given in (5.5) and (5.6). Equation (6.7) is unchanged by the substitution 
(5.5) so that it will suffice to consider real positive values of u for the range p > pz . 


} 4 ‘ ” 1/4 
o> Ba. 3 = KE (p — »2)| > 0. (6.10) 
In contrast to the cases n = 0, 1, Eq. (6.7) does possess solutions of the form (5.6), 


implying that values of p < p.. are possible for n > 2. Equation (6.7) remains the same 
for all u of form (5.6) so that it is sufficient to consider the substitution 


D< Peo; p=t a, n > 0. (6.11) 
Introducing Kelvin functions of order n by the relations 
Jen) = ber,n + tbei,n } 
- L (6.12) 
I,(¢’*n) = exp { bers — ibei,») | 
we write Eq. (6.7) in the form 
; a 
n U,(n) 
ie = i, le 
PSp ef Tn) + (6.13) 
where 
U,(n) = 201 — v)n(n® — 1)§ 2°”’nn(ber, nbei,.n — ber, nbeisn) 
\ 
+ 2'7(1 — v)n(n + |" (bern + beizn) (6.14) 
n 
~= (ber?n + beizn)’ + (ber{n)? + coeis* |} 
n 
and 
T,(n) = [0 — v)’n’(n? — 1) — ']2'*n(ber{nbei,n — ber, nbei;n) 
f | P (6.15) 
+2”*(1 - pat{re| 4 (bern + beitn | — (ber{n)’? — woeitn}. 
| 
Equations (6.7) and (6.13) assume the same limiting form when p = pz , that is 
when yp and 7, respectively, are zero, namely 
Dp = Poa ; “= (1 —v3B4+yn'(r’ - pf + 1(1 — »)(n — 2) 
(6.16) 





n*(n — 1)(1 — »)(4n +9 — ay 
16(n + 2)*(n + 3) ‘ 
7. Solution of frequency equations for n = 2. As the character of the frequency 
calculations is the same for all cases n > 2 we limit ourselves here to the case n = 2. 
According to (6.10) and (6.11) we have that the frequency p is given in the form 


D 
> ° 2 = 3 _— ‘ a 
Dp > Po ; p Po + shat LM, (7.1) 
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D 
SMe 3 p Den = ie (7.2 
Ee I I pha* ] 4.4) 
The quantities yu and » are determined as functions of k* = 48 (1 — v’) (HH /h F by solving 


Eas. (6.7) and (6.13). We carry out this solution by finding «* as a function of u and 7. 
Considering first Eq. (6.7) we find that «° has a value slightly above 25 when u 0 


which decreases to zero when yu be 2.30. For larger values of yw, from u = p,» to 


u = wi, we have «* negative which means that this part of the curve has no physical 
significance. From u us, tou ta, the values of «* are again positive. The next 


branch of the curve fo! positive x lies between the values a and po... The pattern 


repeats itself with increasing values of u (see Fig. 1). 
We note that the values M2.0 > H2.1; Hh , ete., are roots of Eq. (6.7) for x 0, that 


is roots of the corresponding flat-plate frequency equation 


R.(u) = 0 (7.3) 

The values ui , , uf. , ete., are roots of (6.7) for x* ©, that is roots of the equation 
S.(u) — R.(u) = 0. (7.4) 
We note that for k 1, 2, --+ the values of u., and of uw; , are so nearly the same that 


for the purpose of frequency calculations we may disregard their difference. This means 
| | : : 
that for n = 2 and k > 1 (and by implication for n > 2) we may take for practical 
purposes u from the flat-plate frequency equation, just as for the cases n = O and n = 1. 
However, in addition to the regime in which this is possible we now have the regime 
k = 0 where this is not possible. Numerical values of y»; , are given in Table 1. 
My 
We now turn to Eq. (6.13) and determine 7 as a function of «*. We find a curve which 


; P ae ee — é' 
gives x somewhat greater than 25 for 7 0, with increasing values of x° as n increases. 


It is apparent that the n(«*)-curve represents a continuation of the u(*)-curve for k = 0. 
The common point of these two branches for u n = Ois, as it should be, given by Kq. 
(6.16 

On the basis of the data contained in Fig. 1 we have calculated the frequency curves 
n = 2,k = 0, 1, 2 as contained in Figs. 2 and 3. It is apparent that there is an essential 


= 0 where we have no nodal circles and the cases k > 1 


difference between the case k 
0 the influence of shell curvature 


where we have one or more nodal circles. For k 
remains insignificant, the frequency being very nearly equal to the corresponding flat- 


plate frequency. For k > 1 the influence of shell curvature is qualitatively significant. 


TABLE l. Ln tforvy = 1/3. 








k 
| 0 1 2 3 
0 , 9 076 38.52 87.82 
1 — 20.52 59.86 119.0 
2 5.251 35.24 84.38 153.34 
3 | 12.23 52.91 112.0 190.70 
4 21.49 73.9 142.5 231.03 
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8. Frequencies of membrane vibrations. The differential equations of membrane 
theory follow from (2.1) to (2.5) by setting in them 


D=0. (8.1) 


If this is done the differential equations (2.1) and (2.2) may be reduced to the following 


form 
“2 . Cc! 
v"| ok = a. R? w| = Q, (8.2) 
1 Vy = —<» aw. (8 3) 
R ~ at’ 7 


f ) 
; , cos n6@ | ints ; ; 
Equation (8.2), with w = J \ w,(r)e'”', is satisfied either when 


sin nO} 
C 
or, when 
w, = A,r” + Br", l1<n, (8.5) 
Wo = Ay + B, logr. (8.6) 


In either case F’,, is found in terms of w, by integrating the relation 
LF, = phRp’w, . (8.7) 


Since D = 0 only the first two of the four boundary conditions (4.1) and (4.2) need to 
be considered. Appropriate substitution of (8.5) to (8.7) shows that these equations 
are satisfied only when p = 0. 

There remains the frequency p = p. which is associated with arbitrary functions 
w, , insofar as the differential equation (8.2) is concerned. We have already seen in Sec. 
5 that the frequency p. is the limiting solution for all modes n = 0, 1, and k > 1, as 
D — 0, or equivalently as H/h — o. For n > 2 inspection of the solutions (4.7) and 
(4.8) and of Eqs. (6.1) to (6.4) reveals that the conditions 


D-—O, x > @ (8.8) 


are equivalent. Furthermore, for any given mode k > 1, yw as given by (6.7) remains 
finite as x‘ — . Equation (7.1) indicates, therefore, that for modes with one or more 


nodal circles, 


limp — Po ; k > 1. (8.9) 
D-0 


However, since the one root of (6.13) tends to infinity at the same time that x‘ > © 
we have that the frequency p,.o given by (7.2) associated with no nodal circles does not 
tend to the frequency p.. In the above sense we consider the frequency p. as the mem- 
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brane vibration frequency of the shallow spherical shell. Figure 4 which represents values 
of p./Pn.. a8 functions of H/h illustrates the way in which the membrane frequency 
is approached. 

















Fic. 4. Relation between shell frequency p,., and membrane frequency po, when vy = 1/3. 


9. Frequencies of inextensional vibrations. The differential equations for inexten- 
sional vibrations are obtained from (2.1) and (2.2) by setting in these equations 


I 


gy (9.1) 
"& 
It has been shown earlier [3] that inextensional vibrations of shallow spherical shells 
with free edge can occur only for the modes k = 0, n > 2, with frequency p, given by 
H ( er 
= = 0: = | Bai (9.2) 
h P pha*, sin 
H Se be ee ' 
>0; 2= :) [ACL — (ne? — 1yn*}'”. (9.3) 
} pha’, 
In order to see to what extent these results are contained as limiting cases in our present 
calculation we note that 1/C = 0 is equivalent, just as D = 0 is, to 
K = @ (9.4) 


provided H is different from zero. 


Since when x* — © we have p — p. for k = 1, 2, --- it remains to consider what 
happens when k = 0. According to Fig. 1 the frequency equation in this case for suf- 


ficiently large values of «* is Eq. (6.13). As «* tends to infinity 7 also tends to infinity. 
We find by means of the asymptotic expansions for Kelvin functions that (6.13) assumes 
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the following limiting form for large n, 


K~ ni + 41 — vn? — Dn. (9.5) 
Substitution of (9.5) into (7.2) reduces (7.2) to the following form 
D : ' 
op =H — ve — In (9.6) 
pha 
which shows that when / 0, H # 0 then when h — 0 we have p — p, as given by 
(9.3). 
On the other hand, for k 0, h € 0, letting H — 0 in (7.1) results in p — p, as 
given by (9.2). 
Che foregoing considerations indicate that our frequency formula, for k = 0, contains 


as limiting cases both the appropriate flat-plate frequency formula and the frequency 
formula for inextensional vibrations of shallow shells. The way in which the transition 
from one of these formulas to the other takes place is clearly seen by means of the curve 
which in Figs. 2 and 3 is labeled n = 2, k (0. 

10. The energy of stretching and the energy of bending. Qualitatively, in vibrations 
which are approximately membrane vibrations the strain energy of stretching will be 
large compared with the strain energy of bending. On the other hand, for vibrations 
which are approximately inextensional the strain energy of bending will be large com- 
pared with the strain energy of stretching. 

We consider in this section the quantitative aspects of the foregoing statement by 


calculating the ratio 


y 
Vet V, 


’ 


where V, is the strain energy of stretching and V, is the strain energy of bending. For 
shallow shells these quantities are given by the following expressions: 


a oo) oe | ee 2 |! )) 
‘ 7 = Ii") * 9 , -— 
V, | | yA | vit + al + | (2 (t 00 (10 1) 


) |}. dé dr, 


4 


r= fC Blow +20 ~9f(2 2) 
‘agi , (10.2) 


a’w (1 dw 1 ow 
Or (7 or + r° xt) | one. 


We omit the details of the lengthy calculations involved in the evaluation of (10.1) 
and (10.2). The results of these calculations are incorporated in Fig. 5 which shows how, 
for k > 1, the transition from plate to shell is associated with a transfer of strain energy 
from bending to stretching and how, for k = 0, we start with bending energy, encounter 
a certain amount of stretching energy for moderately curved shells, which, with in- 


_ PF (4 OF , 10% 


PF 
or’ \r or r 300° 


to 


creasing curvature, goes over into bending energy. 
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Fic. 5. Ratio of stretching component of strain energy to sum of stretching component and bending 


= 1/3. 


component as function of n, k and H/h, when v 
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THERMO-ELASTIC SIMILARITY LAWS* 


BY 
A. E. GREEN**, J. R. M. RADOK AND R. S. RIVLIN 
Brown University 


1. Introduction. In a book edited by Hetényi [1], Mindlin and Salvadori have 
discussed certain similarity laws for thermo-elastic problems. They were, however, 
primarily concerned with the replacement of thermo-elastic problems by purely elastic 
problems involving dislocations. In the present paper we are concerned with the con- 
ditions which must be satisfied by the physical parameters describing the properties 
of a body and a scale model (scaled both in linear dimensions and .time) in order that 
the stress in the original body may be determined from that at the corresponding position 
and time in the scale model by means of the multiplier M (= rigidity modulus of body/ 
rigidity modulus of model). It will be seen that the conditions for the stresses in the 
body to be an arbitrary constant C times those in the model are then easily obtained. 
In deriving these conditions much of the mathematical analysis of Mindlin and Salvadori 
could have been used. However, we have preferred to obtain them by somewhat different 
methods. 

We consider that surface tractions and temperatures are specified on the surface 
of a body and a scale model, the surface tractions and temperatures on the body being 
respectively 7 and ©, (a constant) times those at corresponding points and times in 
the scale model. The dimensions of the body are assumed to be / times those of the 
model and any time for the body is 7» times the corresponding time for the model. 

It is then found that for three-dimensional thermo-elastic problems, the stresses in 
the original body are MW times those in the model provided that 

i) Poisson’s ratio o is the same for the model and the original body; 
(ii) 8 (= OQ v/u) is the same for the model and original body, where » is the co- 


efficient of thermal expansion, u is the rigidity modulus and 0, = 1 for the 
mode! and 

iii) a (= rox/l°) is the same for the model and original body, where «x is the thermal 
diffusivity and 7» 1,7 = 1 for the model. 


If the body is subjected to plane thermo-elastic strain, the stresses in the body are 
7 times those in the model provided that 8, ¢ and a have the same values for the model 


and the original body. If the surface tractions applied to each closed boundary of the 
body considered are self-equilibrating, the stresses in the body are M times those in the 
scale model provided that 8(1 — 20)/(1 — o) and a have the same values for the model 


and the origina! body. 

In a generalized plane stress problem, the average stresses are M/ times those in the 
model provided that 8, « and a have the same values for the body and the model. If the 
surface tractions applied to each closed boundary of the body are self-equilibrating, 
the conditions become that a and 8(1 — 2c) have the same values in the model and the 
original body. 

*Received November 23, 1956. The results contained in this paper were obtained in the course of 
work on a project sponsored by the Armstrong Cork Co. The authors are indebted to the Armstrong 
Cork Co. for permission to publish the work. 

**(On leave during 1955-56 from King’s College, University of Durham, Newcastle-on-Tyne, England. 
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\. 


Finally, we have the case of a body in plane thermo-elastic strain and a scale model 
in generalized plane thermo-elastic stress, the surface temperatures in the body being 
@, times the average surface temperatures in the model, while the surface tractions in 
the body are M times the average surface tractions in the model, the averages being 
taken over the thickness of the generalized plane stress model. Then, the stress in the 
body is ./ times the average stress at the corresponding position and time in the model, 
provided that o, 8 and a have the same 
However, if the surface tractions applied to each closed boundary of the body are self- 


equilibrating, these cor 


values in the model and the original body. 


ditions can be replaced by the weaker conditions 
a and B*(1 Zo” Bil — 2a)/(1 a), 


where o*, 8* and a* are the values of o, 8 and a@ respectively for the model. 


In all the cases considered, it is assumed that no body forces are applied and that, 
although the surface tractions and temperatures may vary with time, such changes 


are sufficiently slow for inertial forces to be neglected. In the plane problems, the simi- 


larity laws apply, of course, only to the planar stress components. 
In each of the cases discussed, the modelling conditions given above a 
in the body considered are .W/ times those in the model. Since the 


we Can readily deri e con- 


re those tor 


which the stresses 


governing eq ations of th problems considered are linear, 


ditions under which the stresses in the original body shall be any constant C times 


those in the model. For example, let us consider the case in which the original body 
is in plane strain and the model is in generalized plane stress, and the surface tractions 
applied to each closed boundary of the body are seltf-equilibrating. The condition 
B*(1 — 20” B(1 ’¢)/(1 — o) implies that 


0 
N l — 2a 
where N is the ratio of the thermal expansion coefficients for the original body and the 
model. If we now multiply both the surface tractions and temperatures in the scale 
model by .//C the stresses in the original body will be C times the average stresses 
in the model and ©, will be changed to 


C (1 — 2e*)(1 — o 


O05 : 
N 1 — 2 


Similar results may readily be obtained in the other cases considered. 

2. Fundamental equations. We use Cartesian tensor notation and denote a rectan- 
gular coordinate system by y; , where Latin indices take the values 1, 2, 3. Components 
of the displacement, strain, stress and surface traction (or stress vector) are denoted 
by v; , €; , 7; and F; respectively and r is time and 6 temperature. If motion in the 
elastic body is such that inertial terms may be neglected, and if body forces are zero, 
then 

OT 
= (). (2.1) 


The stress-strain relations [2] are 


‘ ‘ Co es " 
T Qule;; + — €,,0;;} — v06;; 
l 20 
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1 / Ov; Ov; 
= HE ™ #1) (2.2) 


and the equation of heat conduction is 


5 a (2.3) 


dy; dy, «x Ot 


In (2.2) 6;, is a Kronecker delta, u and o are the rigidity modulus and Poisson’s ratio 
respectively, and v is a coefficient of thermal expansion. In (2.3) « is the thermal dif- 
fusivity. 

We now express (2.1), (2.2) and (2.3) in non-dimensional form by using the sub- 


stitutions 


II 


7 * lx; ] ; ies lu; ] T 


Tol, €5; = Ci; ’ (2 4) 
tr, =wth,, Fr = pf. , 6 = 41’, 


where | is a standard length, 7) a standard time and ©, a standard temperature. Also, 
x, ,u,,t,e:;,t:;, f; and T are non-dimensional coordinates, displacements, time, strains, 


stresses, surface tractions (or stress vectors) and temperature respectively. Using (2.4), 


Eqs. (2.1), (2.2) and (2.3) become 


i = 0, . (2.5) 


OX; 


—~9 2 :) - 
ti; al. + l Ia €,, 93; BT 5;; , 


(2.6) 
1 (ou, Ou; 
~~ 2 (2 + au) 
and 
aT 1 oT : 
ax, ox; adat (2.7) 
respectively, where 
ok 0 
a= * he and B= Oov | (2.8) 
l Le 


The stress components ¢;; may be eliminated from (2.6) and (2.7) to give the three 
differential equations 
07 u; l 07u; oT 
ae of a nn — §—— a 6), (2.9) 
Ox; OX 1 — 2o Ox; 02; Ox; 

In addition to the above fundamental equations, we have boundary conditions in 
which temperature or temperature gradients, displacements or surface tractions—or a 
mixture of such conditions—are prescribed. For clarity, we restrict our attention to 
boundary conditions which involve no new parameters, although it would not be difficult 


to include these. 
If F, is the surface traction per unit area acting on a surface, the normal to which 
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has direction-cosines 1, , then 
F Fishy (2.10) 


We see that the non-dimensional surface traction f; is given by 


f, = ,l,. (2.11) 


t 


If temperature or temperature gradients 


3. Similarity laws in three dimensions. 
that 


are prescribed on the boundaries of the body considered, then we see from (2.7 
the non-dimensional temperature distribution throughout the body depends only on 


the non-dimensional coordinates x; , non-dimensional time ¢, and the non-dimensional 


parameter a. If we regard this non-dimensional temperature distribution as being the 
actual temperature distribution in a scale model of the body under consideration, then 
the temperature distribution in the original body may be obtained by using the scale 
relations (2.4), provided a is kept constant. If, however, we are concerned only with a 
steady state temperature distribution, we see from (2.7) that 7’ is then independent of 
a, so that no restriction on a is necessary. Further, if the boundary conditions on the 
temperature do not involve a time scale (e.g. if their dependence on time has the form 
of a step-function), then we can choose the scale-factor 7, in such a way that a is kept 
constant. 

If boundary conditions for the body under consideration are given in terms of dis- 
placements, then (2.9) shows that the non-dimensional displacement u; throughout 
the body depends, in general, on ¢ and 8. Again, if the boundary conditions are given in 
terms of surface tractions (i.e. in terms of stress components), we see, from (2.5), (2.6) 
and (2.11) that the non-dimensional displacement components u; and stress components 
t;; throughout the body depend, in general, on o and 8. Thus, if we regard these non- 
dimensional displacement and stress distributions as being the actual displacement and 
stress distributions in a scale model of the body under consideration, then the displace- 
ment and stress distributions in the original body can, in general, be calculated from 
those in the model, by means of Eqs. (2.4), only if ¢ and 6 have the same values in the 
model and the original body. 

4. Plane strain. The assumptions of plane strain are that 


u, = U(X, , 22,8), Uo = UAT; , 2 ;, 8); ot TO 4 Bed; u = 0. (4.1) 
It follows from (2.5), (2.6) and (2.7) that 
€31 = Coz = € = t,, = ft, 0, (4.2) 
and 
Ony = 0, 1.3) 
On 


- (4.4) 
Pa I] [= Ou,, 
= 2 \Or, T OX) ), 
eT laf in 


OX, OX) a al 
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where Greek subscripts take values 1, 2. With the help of complex variable notations, 
Eqs. (4.3) and (4.4) can be integrated in terms of complex potentials (see [3] and [4]). 
We use the notation 


x; + 12 9 D = Uy a. Ws ’ 8 = tis a too ’ ® = ty, _ too 4 iti. ’ (4.6) 


“ 


so that Eqs. (4.3) and (4.4) are replaced by 


i> 30 : 
ae t a = % = 
2 (aD 2D) 
) = & SS — aon —_ [ _ 
et mai + og) — 78F aii 
2 42D 
7 oz’ 


where a bar placed over a quantity denotes the complex conjugate of that quantity. 
From (4.7), it can be shown that 6 and # must be expressible in the form 


Fi 
@= —4 — 9 
Oz 02’ OZ 02’ (4.9) 





Q=4 


where ¢ is a real function (Airy’s stress function). The second equations in (4.8) and 


(4.9) can be integrated to give 
D+ 2 = 401 - of, (4.10) 


where f(z) is an arbitrary function* of z. Using (4.10) and its complex conjugate we 
may eliminate D and @ from the first equations in (4.8) and (4.9) to obtain 


#e | sta ey — M2: 
a f@at+sfe® A — 0) ve (4.11) 
Since ¢ is a real function, two integrations of (4.11) then give 
wae iad B(1 — 2 ) 
g = 2f@ +22 + g@) + G® — —— &, (4.12) 


where g(z) is a second arbitrary function of z and where R(z, Z, ¢) is a real particular 


integral of the equation 





aR 
4 —— 7. (4.13) 


From (4.9), (4.10) and (4.12), we now obtain stress and displacement components in 


the form 


a =4"@ + 7@) - = 7, 
< ee 4B(1 — 2c) dR 
& = —4e7'"@ + g'"@] + BO 22 SE (4.14) 


— (2 — Ae)tle) — 272) — arts) + AL § 200 OR. 
D = (3 to) f(z) zj (2) g (2) + 1 = 02 


* f(z) and g(z), defined in (4.12), also depend on t, but this is not shown explicitly. 





































386 A. E. GREEN, J. R. M. RADOK AND R. S. RIVLIN [Vol. XV, No. 4 


If X, are the components of non-dimensional resultant force (per unit length of 
x; axis) due to the stresses acting across an are AB in the material, in a certain sense, 
then 


»B »B 


X, = fids = | (ti. dz, — ty dx) 
‘ 


si - (4.15) 
~B »B 
Ae = | fds = | (too dx, — ty dx), 
JA JA 
where ds is an element of length of the are AB. Hence, using (4.6), 
»B 
X, + 1X, = | (0 dz — @ dz), (4.16) 
7A 
and with (4.9), this becomes 
oe B 
X, + 1X, = 2i| 22 , (4.17) 
Oz A 


where [ ]; denotes the change of value of the argument as we move along the are from 


A to B. From (4.17) and (4.12), we have 


— 2c) ak |" 
= aR | (4.18) 


X,+ 1X, = 2i| f(0 + zf(2Z) + 92 — — 
; . , l—o OZ |, 

5. Single-valued stresses and displacements. We now restrict attention to problems 

in which the stress and displacement components and the temperature are single-valued 
and without singularities at every internal point of the body. We also assume at present 
that the real particular integral R of (4.13) can be chosen to be single-valued, together 
with its derivatives up to and including the fourth order. This assumption will be justified 
in the Appendix (Sec. 9). It follows from (4.14) that the contributions of FR to the stress 
and displacement components are single-valued. The conditions for the stress and 


displacement components to be single-valued therefore reduce to 


I 


roatfre’w=0, k/’@®+ 9’'@hk = 0 


| 


a 


and 


[((3 — 4o0)f(z) — 2f 


nN 


)— 9D, = O, (5.1) 
where [ ]4 denotes the change in value of the function inside the brackets on passing 
once around a contour in the conventional sense which keeps the area enclosed on the 
left, the contour lying inside the body. 
We may differentiate the last of Eqs. (5.1) with respect to z( or 2), inside the bracket, 

since we have already assumed that the relevant derivatives exist (see [5]), to obtain 

(3 — 40) f’(2) — fala = 0. (5.2) 
Using this, together with (5.1), we obtain 

aa 4 sre A = 

(f’@ls = 0, [g’"(2)Ja = 0, 


(3 — 40) [fla = [HOM , (5.3) 





as given, e.g. by [2] or [3], for the case when the temperature is constant throughout 


the body. 
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If we further restrict attention to problems in which the stress system is such that 
the resultant force acting on any closed curve in the body (or its boundary, if singular 
points on the boundary are excluded by small indentations of the contour) is zero, 
then, from (4.18) we see that the conditions (5.3) may be replaced by 


(f(2]s = 0, [92K =0., (5.4) 


6. Similarity laws for plane strain. Suppose that the body under consideration is 
in a state of plane thermo-elastic strain and that the region S occupied by the cross- 
section of the body is finite and multiply-connected and bounded by one or more smooth 
non-intersecting contours L, , L, , --- , L, , of which L, contains all the others. 

In view of (5.4), we shall assume that, for the particular initial and boundary con- 
ditions of the problem, f(z) and g’(z) are regular functions of z, for every value of the 
time ¢, in the open region S, and continuous in S and on its boundaries, except possibly 
at points on the boundaries at which isolated forces or couples act and where the singu- 
larities are prescribed. We shall also assume that derivatives with respect to time exist 
up to any required order. Then, from (4.14) and (4.18) it follows that the non-dimensional 
stress and displacement components are single-valued and continuous throughout the 


interior of S for all time and that the resultant force acting on any closed contour in S 
is zero. Moreover, if temperatures and self-equilibrating surface tractions are prescribed 
on the boundaries so that the right-hand side of (4.18) is given on Ly , L,, -:- , L, , the 
solution (4.14) for @ and & of the resulting boundary-value problem involves only the 
physical parameters a and 6 (1 — 2¢)/(1 — a). It is thus seen that if we regard the 


non-dimensional stress components given by (4.14) as the actual stress components in 
a scale model of the body under consideration, then the stress components in the original 
body can be calculated from those in the scale model by usingsthe scale relations (2.4), 
provided that a and B(1 — 20)/(1 — oc) have the same values in the model and the 
original body. However, the non-dimensional displacement components D given by 
1.14) involve the physical parameters a, ¢ and 6 (1 — 2c)/(1 — o) and therefore the 
displacement components in the original body can, in general, be calculated from those 
of the scale model only if a, ¢ and 8 have the same values in the model and the original 
body. If the surface tractions on each of the contours L, , L, , --: , L, are not self- 
equilibrating, then the conditions for the stress and displacement components to be 
single-valued are given by Eqs. (5.3). Thus, if the surface tractions on Lp , L,,---,L,, 
given by (4.18), are specified, it is seen that Eqs. (4.14) for 6 and # and Eq. (5.3) involve 
the physical parameters a, B(1 — 20)/(1 — o) and oa. Hence, in order to calculate the 
stress components of a body from those of a scale model, by means of Eqs. (2.4), a, 
o and 6 must, in general, have the same values in the body and the model. 
7. Generalized plane stress. We consider a body bounded by plane faces 7; = + h, 
where h is a non-dimensional constant, and by cylindrical surfaces perpendicular to 
to these surfaces. The planes x1, = + h are assumed to be free from applied stress and 


to be insulated so that there is no loss of heat from them. Thus, 
oT 
OX; 


0, ts: = tos = ty5 = 0 (2X3 = +h). (7.1) 


We assume that the material of the body has Poisson’s ratio o*, coefficient of thermal 
expansion v*, rigidity modulus u* and thermal diffusivity x* and a* and §* are defined 
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.* * 

Tok Oov 

a 2 —_.. de 
at = 2 and p* = ue? 


where / is again a standard length, 7) a standard time and ©, a standard temperature 
and y* is taken as the scale factor for the stress components. 
From (2.6), replacing o and 6 by o* and §*, we have 








* 1 — 20* . 
— = 2 — (t,, + B*T) (7.2) 


= | ar < a 
1 — o* ' 2(1 — o*) 


and using this to eliminate e,, from the remaining equations in (2.6) we find that 
& ©&q 


¢ o* 
hh, = ale, + er. cub | 





_ BL = 2o)Tiry | o*basdry 
1 — o* 1 — o*’ 
ths = 2e3 . (7.3) 


We now assume that the body is subject to a stress system which is such that the 
Uz are even functions of x, and wu; is an odd function of x, . We then 


displacements 4, , 
using con- 


follow the usual procedure and take average values of Eqs. (2.5) and (7.3), 
7.1). We add the extra assumption that the average value of ¢,; is neglected 


ditions 
Denoting the average values of f,, , ¢,, , 


in comparison with the average values of th, . 
u, and T by t,* , e,* , u* and 7* respectively, we obtain 


Ot“. ~ 
ar. = 0, \ , t) 
Ou 
* * ) 
o re B { l —_— wae 
i,* 2 - + _ $5.) Sh ey eee T* by, ’ 
l—o l o 
_ on n ou’) (75 
€ —_ ° o fo) 
2 OX, OX 
and 
oT* 1 o7* "6 
een a eer meets (4.0 
OX, OX) a™ dal 


If we denote by X* the average values of the components of the non-dimensional 
resultant force (per unit of z-axis) due to the stresses acting across an arc AB of the 
material lying in the 2z,z,-plane, then we see from (4.15) that 
P »B 
Xt= | (tHdx,—ttdz,), X=] (ttde, — thr). 


JA 


replacing » by u* in Eqs. (2.4), that the average values of the actual 


We note, 
stresses 7*, , displacements vt and temperature @* can be obtained from the non-di- 
mensional values by means of the relations 

Yy lx. , >) ae T= 7), 7). = y*t* . 


and 


% 
rl 
® 
* 
~I 
90 
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Equations (7.4) to (7.7) are identical in form with Eqs. (4.3), (4.4), (4.5) and (4.15) 


if we replace 


te = Ox T*, wut, oF/(1—o*), B*(1 — 20*)/(1 — o*) 


? 
and a* in the former by 4, , & , T, uw, ¢/(1 — 20), 8 and a@ respectively. Hence, we 
can write down the solution of Eqs. (7.4) to (7.7) by analogy with the plane strain 
problem. The conditions for the average values of the non-dimensional displacement 
and stress components to be single-valued may similarly be written down by analogy 
with the case of plane strain. 
With the notation 
ee = (i+ tt, ~* = 14 — t$+ 2t, and D* = uf + it, (7.9) 


we obtain [see Eqs. (4.14)]} 
O* = 4[f*/(2) + f*’®] — B*(1 — 20*)T*, 











@* = —4[zf*’"(z) + g*’"(2)] + 46*(1 — 20*) oR* 
: ; 02 02’ 
~ a=, on iaiibd ok* — 
D* = Ido f*(z) — zf*'(2) — g*'(2) + B*(1 — Qo* rae (7.10) 
where f*(z) and g*(z) are functions of z which depend on the boundary conditions and 
R*(z, Z, t) is a real particular integral of the equation 
a°R* , . 
4—— = T*. (7.11) 
Oz OZ 


From (7.7) and (7.10) it is seen that 





02 A 


. , ak* |? 
X% + iX$ = 21} f*(z) + zf*"(2) + G*'@ — B*(1 — 20%) (7.12) 


It can readily be seen, in a manner similar to that adopted in Sec. 5, that the con- 
ditions for the average values of the non-dimensional average stress and average dis- 
placement components ¢,* and u* to be single-valued are 


oO =0, [g*t’"@lk =0 
and 
3 — o* ree ‘ = 1x) 74 4 
1d o* [f (2) ]‘4 = (9 (2) |4 ? (7.13) 


In the case when the stress system is such that the resultant average force acting on any 
closed curve in the body is zero, the conditions (7.13) become 


[f*(2)]4 = 0, (g*’(2)]4 = 0. (7.14) 


8. Similarity laws for generalized plane stress. We now assume, as in the case of plane 
strain previously discussed, that the cross-section of the body considered, normal to 
the x,-axis, is finite and multiply-connected and bounded by one or more smooth non- 
intersecting contours Ly , L, , --- , L, , of which Z, contains all the others, and that 
the values of the non-dimensional average surface tractions (and hence of X*) and 
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non-dimensional average temperature 7* are given over each of the contours Ly , L, , 

L. . Since, from Eqs. (7.6) and (7.11), the values of 7* and R* throughout the 
body depend only on the physical parameter a*, the expressions (7.10) and (7.12) for 
@*, &* and X* involve only the physical parameters a* and 6*(1 — 20*). Also, if the 
surface tractions applied to each of the contours Ly , L; , --- , L, are self-equilibrating, 
the conditions (7.14) for the average stress and average displacement components to 
be single-valued do not involve any physical constants. It follows that, in the case when 


the specified surface tractions applied to each ol the contours Ly , L, , +: | DL, are self- 


equilibrating, if the non-dimensional average stress distribution is taken as the actual 
average stress distribution in a scale model of the body under consideration, then the 


average stress distribution in the original body can, in general, be calculated from that 


in the model by employing the relations (7.8), 11 a” and 8*(1 — 2c*) have the same values 
for the mo le] and the o iginal body. In the case when the specified surface tractions 
applied to each of the contours L, , L, , +: , L, are not self-equilibrating, since the 
conditions (7.13) for the stress and displacement components to be single-valued involve 
o*. the average stress distribution in the original body can, in general, be obtained 
from that in the model by means of the relations (7.8) only if a*, o* and B* have the 
same values in the mode! and the original body 


Now, suppose that we have two bodies A and B with geometrically similar cross- 
sections bounded by the « osed contours La , y = eae L, . as de scribed above. Let the 
dimensions of the cross-section of A be / times those of B. The body A is maintained in 
a state of plane thermo-elastic strain and B is maintained in a state of generalized plane 
thermo-elasti ‘ stress bv surface tractions and surface temperatures such that the actual 


nd average surface temperatures in B are equal to the non- 


J 


average surface 
dimensional surface tractions and surface temperatures in A. We assume that the 
physical parameters for the body A are the unstarred parameters defined in Sec. 2 
while those for the body B are the starred parameters defined in the present section. 
Taking / llr Lp 1 and @, = 1 for the body B, we see that (7.10) and (7.12) 
give the actual average stress components and surface tractions in B, while 7* represents 
the actual average temperature in B. Then, comparing Eqs. (7.6), (7.10), (7.11), (7.12) 
and (7.13) with Eqs. (4.5), (4.14), (4.13), (4.18) and (5.3), and bearing in mind that 
we have assumed 7 T* and X , = X* on the surfaces of A and B, we see that provided 


that 


ies 
| 
Ww 
Q 
| 


at =a, B*(1 — 2e*) = ——— : = 3 — do, (8.1) 
l—o 1 + o* 
i.e. provided that 
fe Co B(1 — 2a) 
ra’ a c* = and £* ~ ae (8.2) 
l-—o l — 3e 
we have 
6*=060 and &* = 9, (8.3) 


so that the stress distribution in A can be calculated from the average stress distribution 
in B by the scale factors defined in Eqs. (2.4). 

If the surface tractions applied to each of the contours Ly , L, , --- , L, are self- 
equilibrating, then comparing Eqs. (7.6), (7.10), (7.11), (7.12) and (7.14) with Eqs. 
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(4.5), (4.14), (4.13), (4.18) and (5.4), we see that Eqs. (8.3) are satisfied provided that 


BU — 20). 


l-—o 


a* =a and £6*(1 — 2c*) = (8.4) 

9. Appendix. The function R. In deriving in Sec. 5 the conditions for the stress and 
displacement components to be single-valued the assumption was made that the function 
R(z, 2, t) satisfying Eq. (4.13) can be chosen to be real and single-valued, together with 
its derivatives up to the fourth order. This assumption will be justified in the present 
section. 

We deal first with the case of steady heat flow and statical equilibrium for which 
ar T(z, 2), 8 = (2, 2) and 

IT» d Ro 


= 0 T= er ov. 
Oz 02 ° ' : Oz 02 (9.1) 


Since 7’, is real, the first of Eqs. (9.1) may be integrated in the form 
T, =W(®+h, (9.2) 
where, since 7, and its derivatives are single-valued inside the region S occupied by the 

cross-section of the body considered, 

[h’'(2)}4 = 0, [h'(2) + WA = 0. (9.3) 
We now suppose that the region S is the multiply-connected region bounded by 
contours Ly, , L,; , --: , L, , as described in Sec. 6. We assume that the temperature 
distribution in the body is such that h’’(z) [which is seen from (9.3) to be single-valued] 


is a regular function of z in the open region S. From Laurent’s theorem it can be seen 
that h’’(z) must be expressible in the form 





- a, b.+ id, i 
h'(z) = ee ce ons of &'"(g), (9.4 
2 2Qn(z — 2%) p> 2n(z — 2z)° ) 
where k(z) is a regular function of z in S, b, and d, are real constants, z, is a point inside 


the contour L, (and therefore outside S) and a, is a constant. The result (9.4) follows 
from the fact that if terms of higher degree than the second in 1/(z — 2) are included 
in the Laurent expansion (9.4) for h’’(z) they can be absorbed into the function k’’(z) 
without vitiating the regularity of k(z) in the region S. 

By integration of (9.4), we obtain 


".. 7 ~ ay Coes Oe = b, 4- id, v Q* 
h'(2) p> 5. log (@ — 2) + > gage + k’(2), (9.5) 
he) = + ( — a) log (2 — 2%) + Sy ett hog (@ — 2) 
— — (9.6) 
oe = ue, _, 
+ k(z) do os (2 — 2). 


From the second of Eqs. (9.3), it is seen that a, is a real constant. It is apparent from 
(9.2) and (9.5) that if h(z) is expressible in the form (9.6) with a, real, 7, is a real, single- 
valued function in S. 
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Introducing the expression (9.2) into the second of Eqs. (9.1), we obtain 


ra) R + Pane ~ 
{| —— h'(z) + h’(2). 9.7) 
OZ OZ 
Introducing the expression (9.5) for h’(z) into (9.7), it is apparent that (9.7) has a par- 
ticular integral for RP, of the form 


; = @ 5 
4R, Da (2 2.)(2 — 2,)[log (e — 2.) + log (2 — 2,)] 
aa (9.8) 
: b 1. ad b. — id . 
F k - F i Z - , = 
-{- p> ! a * {- — s [log (e — z,) + log (2 — 2,)] + U(z,2), 
where /(z, Z) is a real single-valued function in S, given by 
= - se . Cy = \ 
U(z, 2) = zk(z) + 2k(2) — D> 22. (9.9) 


Since the expression on the right-hand side of (9.8) is real and it and its derivatives 
are single-valued, 2, may be chosen in such a way that it is real and it and its derivatives 
are single-valued. 

When the heat flow is not steady, T and RP satisfy the equations 


oT 1 oT oR Ma : 
> = tenor 0, (——— = 7’. (9.10) 
OZ OZ a ol OZ OZ 
Eliminating 7', we obtain 
o / aR 10R 
mea: Nigel —} = 0. (9.11) 
Oz OZ Oz OZ a dl/ 
Equation (9.11) is satisfied by 
oR ] oR q 
4———-—=4Q, (9.12) 
OZ OZ a al 
where @ is any real function which satisfies the equation 
OW < 
——— a (), (9.13) 
OZ OZ 


We make a particular choice of the function Q that it be expressible in the form 


l oP 
gu 12 
a ol 


+7, , (9.14) 


where 7’, is a real single-valued function, independent of /, which satisfies the equation 


ols ie 


Oz OZ 


0 (9.15) 


and P is a real function which satisfies the equation 
oP 


02 02 


= (0. (9.16) 
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Since 7, satisfies (9.15) and is real and single-valued, we see, as we did in discussing the 
function 7, given by the first of Eqs. (9.1), that there exists* a function R)(z, Z) which 
is real and single-valued together with its derivatives up to the fourth order and satisfies 
the equatio1 





T = 4 a Ry (9.17) 
Oz Oz 
Defining R,(z, Z, t) by 
R=P+R,+8,, (9.18) 
we see, from (9.12), (9.14), (9.16), (9.17) and (9.18), that R, must satisfy the equation 
soe 1S ~ (9.19) 


Substituting from (9.18) in the second of Eqs. (9.10) and employing (9.16) and (9.17), 


we have 


T=T,+7,, (9.20) 
where 
a aR, f ; 
mn=4 Oz Oz “— 
From (9.19) and (9.21 
T, = 1 oft, (9.22) 
a dl 


Since 7 and 7, are real and single-valued, we see from (9.20) that 7, must be a real 
single-valued function of z and 2 for every ¢. It follows from (9.22) that we may choose 
R, to be a real single-valued function of z and 2 for all ¢. It must, of course, also be single- 
valued in ¢ and possess derivatives with respect to z, Z and ¢ of any required order. We 
have already seen that R, can be chosen to be real and single-valued and P has been 
chosen to be real and single-valued. It follows from (9.18) that R is real and single- 


valued 
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BOOK REVIEWS 


(Continued from p. 340) 


Handbuch der Laplace-Transformation. By Gustav Doetsch. Volume II. Verlag Birk- 
haiuser, Basel and Stuttgart, 1955. 436 pp. $13.15. 


This book is the second volume of the author’s trilogy on the Laplace transformation. Volume I, 
which appeared in 1950, was concerned with the theory of the Laplace transformation. Volumes IT and 
III are dedicated t6 the applications of the Laplace transformation in various domains of pure and 
applied mathematics. The present volume is divided into three parts, entitled: Asymptotic expansions, 
Convergent expansions, and Ordinary differential equations. An introductory chapter collects together 


the main results and “rules” concerning the one sided Laplace transformation f(s) = |o e~*F(t) dl 
and the two sided Laplace transformation f(s) = } e~*t F(t) dt. The first part centers around the 
following general definition of an asymptotic expansion (“im Sinne von Poinecaré’’): a function g(z) is 
said to have an asymptotic expansion } Br .o Cjh;(z) in a neighborhood U of z provided that for each 
n = 0, 1, one has g doha0 cjhj(z) = O(hn(z)) as z — 2. (It is supposed that ha(z) ¥ 0 in U.) 
Given any functional transformation 7’, (so that 7'(F) = f, say), a theorem which concludes something 


about the asymptotic behavior of the image function f from a knowledge of the asymptotic behavior 
of the original function F is designated as an “‘Abelsche Asymptotik’’. If the réles of the image function 
and the original function are interchanged in the last statement, then such a theorem is known as a 
‘“‘Taubersche Asymptotik”’. The main portion of the first part is dedicated to the Abelsche Asymptotik 
of the one and two sided Laplace transformations, the Mellin transformation, and the complex inversion 
integral of the Laplace transformation (this last is of particular importance in technical applications, 
when the properties of a function defined by an inversion integral have to be determined). There is also 
a chapter on the Taubersche Asymptotik. The second part devoted to convergent expansions, consists 


of two chapters, the first on faculty series 


> a,n! 
eo Se ieee 


s(s + 
) ‘ | 


(s + n) 


and the second on various series expansions, e.g. theta function, Laguerre polynomials, and others. 
The main concern here is with convergent series which are obtained by term by term application of 
the Laplace transformation. The first chapter of the third part deals with the ordinary differential 
equation with constant coefficients (which may be written in symbolic form thus: p(d/dt)Y = F, where 
p(x) = 2" + Cait . + ¢2 + cyand F = F(é)) on the interval 0 < t < ©, the initial conditions 
being 
r. lim Y(),---, Yo" = lim Y""'"'(0, 
na + 

while the next chapter considers the same differential equation on the infinite interval <t< + @ 
with various boundary conditions at — © and + o. There are detailed applications to the theory of 
wave filters and electrical networks, among others. The system of ordinary differential equations with 


constant coefficients 


pi(D)Y; + pil(D)V2 + +++ + pixv(D) Vy = Fb) 


Dyi(D)Y, + py D)Y2 + ++: + pwn(DYx = Feld, 
where 


p(s) = cits” + ++» + ei*s + ci j,k =1,---,N) 


is also treated (specially in the ‘“Normalfall’’). The remaining chapters develop the theory of ordinary 
differential equations with variable coefficients. As is his custom, the author has included a large number 
of valuable literary and historical remarks at the end of the volume. This book is remarkable for the 
wealth of detail and the clarity of the exposition. 

J.B. Diaz 
(Continued on p. 435) 
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THE RECTIFICATION OF NON-GAUSSIAN NOISE* 
BY 
JAMES A. MULLEN (Raytheon Manufacturing Company, Waltham, Massachusetts) 
AND 
DAVID MIDDLETON (Concord, Massachusetts) 


Section I. Introduction. This paper considers the effects of a general ciass of 
non-gaussian random processes in an a-m receiving system. Analysis of the effects 
of non-gaussian noise is needed both for use in problems where the noise is known not 
to be gaussian, e.g. some kinds of radar clutter and atmospheric static, and to indicate 
in uncertain cases how critical the assumption of normal statistics may be. 

The probability densities of a generalized Poisson process and a new approximating 
series for the densities are presented in See. II, with discussion of and results for the 
rectification problem following in See. ITI. 

The asymptotic approximation to the probability distributions, using derivatives 
with respect to the second moments, is new in its general form, although foreshadowed 
by some special results obtained earlier by Edgeworth [1] and Pearson [2]. The problem 
of rectification of non-gaussian noise does not appear to have been attacked with reason- 
able generality before. A number of papers have dealt with the deviation of rectifier 
outputs from normal statistics when the input is normally distributed [3, 4, 5] and some 
results are available for non-gaussian noise in quadratic detectors [6] and for very narrow 
post detector filters [7]. We have obtained expressions for the output covariance function 
of a half wave vth-law rectifier with a sine wave carrier and non-gaussian noise input. 
Explicit results are presented for linear and quadratic detectors. A qualitative discussion 
of the behavior of rectified non-gaussian noise for general values of v is also included. 

In general, the work has shown that the difference between the effects of gaussian 
and non-gaussian noise of the same input power is fairly small unless the signal is weak 
or the noise is strongly non-gaussian. Usually non-gaussian noise produces a greater 
output intensity and a broader output spectrum than gaussian noise with the same input 
power and input spectrum. 

Section II. Probability densities. 2.1 Introduction. The probability concepts and 
notation used here are briefly presented in this section. Probability densities are used 
rather than cumulative distribution functions; in the usual way, for example, one writes 
Wily, . Yo; t,) dy, dy» = the joint probability for a stationary random process that 
y will lie in the interval (y, , y: + dy,) at time ¢, , and in the interval (y , yz + dy) 
at time / 

The distribution can also be characterized by its moments (whenever these exist), 


e.g 


Mmnll) = (yY> [| YiY2Wlyi , Yo 3 0) dy: dye - (1) 


*Received December 27, 1956. This work was supported in part through Cruft Laboratory, Harvard 
University, jointly by the Navy Department (Office of Naval Research), the Signal Corps of the U. 8. 
Army and the U.S. Air Force under ONR Contract N5ori-76, Task Order 1, and was carried out mainly 
during the period July, 1953-January, 1955. 

“~The angle brackets here and subsequently denote the statistical average. 
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The Fourier transform of the probability density is the characteristic function, 


F.(f, , & 5 t) = (exp (26:41 + f242)), (2) 
whose power series expansion generates the moments, 
a (7é "(rk n 

Pu (é oa e co $1 oe S27 (Q 

Fié,@;0= > 2 malt) , (3) 


m,n: 


(again, whenever u,,,,(¢) exist). An equally useful set of parameters are the semi-invariants, 


Aewatt); defined by 


ee (2€,)"(2E2) 
log F.(é, , & ; t) 2. ee Mae (4) 


ay . m!n! 


2.2 The Poisson ensemble. The non-gaussian random processes which we shall 
use are Poisson ensembles, i.e. they are composed of sums of independent variables with 


common distributions and uniformly distributed times of occurrence, t/ , viz., 
> ot — tt). 5) 


Physically, the v; represent pulses produced randomly by a noise source. The above 


V(t; K, {ti} 


assumptions imply that the number of pulses, AK, occurring in a long time interval, 7’, 
is a random variable with a Poisson distribution [8] so that A, as well as the set of U , 
is an ensemble parameter 

The characteristic function of the process, for a finite interval of time 7, 
readily derived by taking advantage of the independence of the pulses and performing 
a weighted sum over a random walk for each fixed K [9, 10]. Thus we write 


can be 


. 7. nT)* , — > : K 
f v\S »f2 5 Ur ; ar ae "[F2 8 5Orl', 


- (6) 


= exp {nT [F.,(& ,& ; Or — 1}, 
where n is the average number of pulses occurring per unit time. F., is the characteristic 
function for the variable v, representing an individual pulse, and F,, is the characteristic 


function for the sum, V. 
The finiteness of the time interval can be eliminated by first expressing the exponent 


as an average, 


oj; be — 1] = nT(exp (tév, + 2&0.) — 1), (7) 


nT [F >. 


itr 
Str 


in which the average must be carried out over all the random variables of a typical 
pulse, e.g. amplitude, shape, phase, duration, or time of occurrence. Let us maintain 
the average sign as a reminder of the presence of other random variables and explicitly 
consider the averages for the times of occurrence. If the time variables are transformed 
from t{,tjtot—u,& ’ then, taking advantage of the uniform distribution of the 
occurrence time, we have 
ae : l t— t 
W(t{ , t) dt{ dt} = W(t — t) 7? a(—4) d(t; — th), (8) 


where 7 is the mean duration of a pulse, so that letting 7 increase indefinitely, we get 
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; . ; t, : 

Foy(é: , 2 5 t) = exp {> / ({exp [2&0 (to) + tEv(t + Y)} — 1) a(')\ , (9) 
—o ¥, 

where y = nr is the average noise “density,” i.e. the average number of pulses per 


second multiplied by the average duration of a pulse. 
The semi-invariants of V are readily obtained by expanding the inner exponential 
ol Kg (9), yielding 


Ann(t) = yv(te)"v(to + 1"). (10) 


One is at once tempted to say, from Eq. (10), that the semi-invariants of V are 
proportional to the moments of v. This is not entirely proper, because the presence of 
the —1 under the average in Eq. (9) acts as a convergence factor, so that the exponent 
is not a true characteristic function, and hence the right hand side of Eq. (10) is not a 
true moment. For a finite time interval, no convergence problems arise, and Eq. (6) 
shows that then indeed the semi-invariants of V are proportional to the moments of v. 

The parameter of the Poisson ensemble which has most influence in determining the 
general features of the noise is the density y. As the noise density increases, the noise 
distributions tend toward the gaussian. When the density is small, the noise pulses 
overlap only slightly so that the noise has the attributes of a deterministic, interfering 
signal, whose time-structure is essentially that of a single typical impulse. These two 
limiting regions are important, since the exact distributions are usually too complex to 
be used analytically, and must be approximated differently in each of the two cases. 

2.3 Narrow band noise. Here, in the course of rectification, the detector or rectifier 
is preceded by a frequency selective network which passes only a band of frequencies 
narrow compared to the center frequency of the band. This is simultaneously desirable 
in order to discriminate against unwanted signals in other frequency bands, and more 
or less inescapable because of the inherent characteristics of the elements of which the 
receiver is built. Thus, although the noise at the input of the receiver may be, and usually 
is, of relatively constant strength over a wide band of frequencies, the noise presented 
to the detector is narrow-band because of its passage through the frequency selective 
parts of the receiver. ; 

Important simplifications of the probability densities can be obtained by taking 
explicit cognizance of the narrow-band nature of the processes with which we are con- 
cerned. A narrow-band random variable can be expressed as 


V(t) = R(2) cos [wot — 0(d)], (11) 


where w, is the central frequency of the noise spectrum and F and @ are random variables 
whose variation with time is slow compared to that of cos wot. Thus, if one calculates 
the moments of V using time averages (over a single member function of the V ensemble, 
which is assumed to be ergodic), R and @ can be assumed to be constant over a single 
period of cos wet. Consequently, all the moments of V which are of odd degree vanish. 
The constancy of envelope and phase over one period of the center frequency is, of course, 
an approximation and cannot be expected to be valid for moments of arbitrarily high 
degree, since a change in R of 6R will give a change in R" of nR** 5R, and thus, for 
sufficiently large n, will be comparable to the change in cos wot. In the approximating 
distributions which will be used, however, only low-degree moments will appear, so 
that the invariance of the slowly varying parts over a single period of cos wot is a valid 


assumption. 
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The presence ol high lrequency terms, such as cos wf, In the expression for the 
second-order moments is not eliminated, by any means. In order to indicate clearly 


the presence of high frequency terms, let us define an envelope factor for the moments, by 


Mns(t) = tim + + 2) opapmy (12) 
P(1/2)({m + n]/2 + 1) si 
Let us further assume that the phase change of the noise voltage, i.e. 6. 6, , is not 
arandom variable*. This is not a valid restriction on some types of random waves; but, 
in this paper, we shall need the results below in their application to individual pulses 
[the of | q. (5)], where the phase would be expected to be constant. 


The second-order moments can now be written 


uy (e R, cos (apt Of, COS (Wolo + wol 10; + al) 
RR. cos” (wot i] COS (Wyl a 
RR cos (wol #,) sm (wyl 6,)) sin (wol a) 
4 R R COS (wal a 
\/ COS (Weyl a (13a) 
and similarly 
us (t V5,(0) cos (at — a (13b) 
Moll M,.(0[} + §% cos” (wot — a)], (13¢) 
which, with y,,,,.(0 u , defines all the non-zero moments up to the sixth degree. 


We note that by comparing the two expansions of the characteristic function, qs. 


(3) and (4), a similar separation of high and low frequency factors can be accomplished 


for the seml-invariant e.g 
Ail A,,(Q) COS (wol a), (14) 
and that a special notation is usually employed for the autovariance function, i.e 
{ R(t yr(t Yro(t COS (Wol a), (15) 
where 7(0 “(0 


2.4 Nearly gaussian distributions. A nearly normal distribution, e.g. Eq. (9) when 4 
is large, can be exp nded asymptotically in terms of the limiting gaussian form. For the 
case where the odd degree semi-invariants are zero; the leading terms of the expansion 


of the characteristic function are 


\ A3,( 0) \o0(t A, a(t ' Aos 4 
i; j ‘ + Ss ey 3 + a {+ 4 4 -7 
s > ' "7 sis >! »r o's 7 $is2 {! Se ' 
t e yA } 
{ ) (16) 
/ 
} ¥ | Orr +2 
‘exp 4 — 5 (6: + 26&r() + £2] ?- 


— 


n taking the Fourier transform of this equation, we note that the multiplications by 


£, will become differentiations with respect to V; , and the two-dimensional form of the 


well-known Edgeworth series is obtained. The terms containing the fourth degree semi- 


*Formulae for the case where the phase change is random may be found in Mullen and Middleton [11]. 
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invariants can be shown to be of order Ay/~* ~ y ' and the neglected terms to be of 
order 7 

{n alternative form, whose derivation is given in Appendix A, can be obtained in 
which derivatives are taken with respect to the second moments, 


Pir. .t:0 1+L+ Oy ’)] 
‘exp [- bl yiki + 2§,&.Yr(t) cos (wot — a) + P.to]} ly. =¥a~¥ 
where 


t Az,(b) 
3! dy 3! OW, Ar 


A, ra)  @ 0 
> 
) 


p Aw? FE rao -| ee | ae 
3! Ors OW, Wo 3! = Ayoro 3! dy, 
and subseripts have been placed on the second moments to enable each term in the 
exponent to be differentiated separately. 

Only semi-invariants of low degree appear, so that the narrow-band approximation 
is still valid. Equation (17) is applicable to any narrow-band, nearly normal distribution, 
whether derived from the Poisson ensemble or not. 

l’or moments after a non-linear operation, the comparative simplicity in the analysis 
when the noise possesses a gaussian distribution is a strong point in favor of the series 
in parametric derivatives. Applying the differential operator to the result for (nonsta- 
tionary) gaussian noise is likely to be involved, but is certainly straightforward. Further- 
more, Where the narrow-band structure of the output is important, this characteristic 
function possesses the additional advantage of presenting the high-frequency part only 
once (in the exponent) while the Edgeworth form has high-frequency terms of different 
orders scattered about among the various semi-invariants, as well. 

2.5 Low-density distributions. When the density is small, the exact form of the 
distribution is too difficult to use, and at the same time, the Edgeworth type of approxi- 
mation is no longer applicable. However, a useful ascending power series in y can be 


obtained, starting from Eq. (9). Let 


; 4 ; ; l 3 
f(& ,&3;0 =1+ | (exp (2&0, + 2&0.) — 1) a(4) , (18) 
so that 
Fifi; yt 30 = exp fh ,&3;0-M =e" HUse,& 50", (19) 
mo fF . 


a form which much resembles Eq. (6). 

The utility of Eq. (19) lies in the fact that the first two terms can be used in their 
exact form and the remaining terms are similar to nearly gaussian distributions, where 
n corresponds to the noise density. After considerable algebra, one obtains 


; , om fo) 7" J y’ a 
y. ft ey 7. £2 - -« a — — 
F,, S1 9» $2 3 t ¢ + ¥é mAs 9S2 3) t) + »» ni‘ | Qn oy” 


yh. (20) 


( 
“exp {_ #y l&} + 2&,E.ro(t) cos ot — a) + 2) + 
2y ; Sais "s n 


-= [Viti + 2WéEro(t) CoS ol — a) + vei} : 
2y leedss6s 


‘exp 4 


——_ 
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where L is the differential operator that appears in Eq. (17). Since our aim is to obtain 
the covariance function after rectification, which in our representation of the rectifier 
(vide Eq. (22) infra), is itself a moment of the input distribution, although in general 
not one of integral degree, an Edgeworth series in which the low-degree moments are 
repre duced exac tly is ppropriate. As a fit to the distribution itself, however, the approxi- 
mation cannot be « X pe ted to be nearly SO good. 

2.6 Example. The first-order probability distribution of rectangular c-w pulses 
will serve to illustrate the behavior of the distributions of the process as the density is 


varied. 

The distributions corresponding to the terms of fixed n in Eq. (19) have been pre- 
viously calculated [12]*. From these, one can easily obtain results, which are shown in 
Fig. 1, for the Poisson ensemble together with a nearly gaussian approximation for 
y = 1 and a gaussian distribution. The figure is a plot of the probability that the value 
of the abscissa, measured in units of the standard deviation, will be exceeded in absolute 
value. 
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Fic. 1. Non-gaussian probability distributions. 


When the noise density is one, the nearly gaussian approximation is a fairly close 
fit to the tail of the distribution, but is significantly different for small values of V. 
With the Poisson ensemble there is a finite probability of having no noise at all, so that 
*Just as in Eq. (10) the right-hand side is not properly a moment, so in Eq. (19), f is not properly 
a characteristic function. However, the appearance of the anomalous behavior is entirely in the time 
structure of the process so that, for the first-order distribution, each term of the power series in y can 


be considered as a characteristic function. 
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the non-gaussian distribution should contain a 6-function of strength e~” at the origin. 
The nearly gaussian expansion, because of its asymptotic character, is unable to repre- 
sent discrete probabilities.* 

When plotted against V, the low-density curves lie above one another; but dividing 
this abscissa by the standard deviation, which varies as y'’*, causes the curves to intersect. 
As the density decreases, the probability that the noise will be zero becomes larger. 
This in turn decreases the standard deviation faster than it does the tail of the distribu- 
tion, so that the probability of relatively high values increases, for smaller noise densities, 
relative to the larger values of noise density. 

Section III. Properties of the rectifier output. 3.1 Introduction. Our next task is 
to use the distributions found in Sec. II to determine the autovariance function of the 
output of an a-m receiver when an additive mixture of c-w signal and non-gaussian 
noise is impressed on the input, viz.; 


Vin() = Vxlt) + Ao COS wol. (21) 


The receiver in an amplitude modulation system contains a band-pass filter which 
eliminates all spectral components except those in a narrow-band around the central 
frequency to which it is tuned, and a demodulator, consisting of a half-wave vth-law 
detector followed by a low-pass filter. The relation between output and input of a non- 
linear device is the dynamic transfer characteristic. For the half-wave (zero memory) 
vth-law device, it is 

IBV’, } > 0 (99 


I= g9(V) = 22) 
0, Vv <0 
where J is the instantaneous output, V the instantaneous input, and @ is an appropriate 
proportionality constant. The non-linear device produces an output whose spectral 
components lie in zones around multiples of the central frequency of the input [14]; the 
ideal) low-pass filter eliminates all these zones except the low-frequency one correspond- 
ing to the zeroth harmonic. 

We wish now to find the autovariance function (or since the input process is assumed 
to be ergodic, the correlation function) of the output. This is by no means the most 
complete statistical description that one could wish for, although the correlation function 
yields considerable insight into the nature of the output. However, to obtain the prob- 
ability distributions of the low-frequency zone is a very difficult problem, which has 
been solved only for the quadratic detector with gaussian input noise (since no moments 
higher than the fourth order are then required) [15, 16]. 

Some general statements about the output covariance function can be made immedi- 
ately. The covariance function equals the mean square of the output when its argument 
is zero. As the argument, ¢, increases beyond bounds, the two functions to be averaged 
become uncorrelated, so that the covariance function becomes the square of the output 
mean. The behavior peculiar to R(t) is revealed by its variations with time, which can 
best. be separated by defining a normalized output covariance function 


R(t) — Ri oo ) 


= : 9° 
~ RO) — R(o) (23) 


Tout) 
*A similar situation arises in the discussion of the Brownian motion, where a 6-function appears in 
the exact form of the velocity distribution which the usual approach of Fokker and Planck is unable 


to provide [13 
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3.2 Noise models. To obtain the output covariance as a specific function of time, 
a definite time dependence must be assigned to the semi-invariants of the noise. An 
immense variety of types is naturally possible; three models felt to be important are 
treated here. In these models the entire time dependence is contained in the envelope of 
the pulses, although in general, variations in amplitude, shape, and phase of the individual 
pulses will also affect the time dependence. Variations in the amplitude or shape appear 
in the semi-invariants in the same way as the variations of pulse envelope, so that no 
loss in generality is to be expected if we attribute the overall time variation to a single 
cause rather than to a mixture of several. This is not true for the phase variation, how- 
ever, since phase changes are associated with the high-frequency part of the semi- 
invariants, and therefore affect the noise distributions in a manner different from the 
amplitudes or pulse envelopes. Most noise sources produce pulses of constant phase, 
so that further refinements are unnecessary; however, the commonest exception, moving 
clutter in a radar system, is certainly an important one.* 

Pertinent data on the three models are summarized in Table 1. The exponential 
pulses represent impulses passed through a single-tuned circuit. Although a single 
selective element is not an accurate model of the tuned stages of a receiver, this type of 
time dependence is important because it is necessary (and sufficient) if the process is to 
be Markoffian in the limit of increasing density. A representative physical case would 
be that of atmospheric static interference in a crystal video receiver. 

The pass-band of an actual receiver is an involved function of the number of stages 
and the coupling networks between stages. As an abstraction from the details associated 
with any particular 7-f strip, we shall take a pass-band of gaussian frequency response 
which preserves the essential features of a pass-band while remaining analytically 
tractable. Admittedly, a gaussian pass-band is not physically realizable, but it is a good 
approximation to the magnitude-frequency curve of actual amplifiers (if not to the 
phase frequency curve), and possesses the important virtue of simplicity**. Impulses 
passing through this 7-f amplifier will become gaussian pulses, as in the second model 
of Table 1. ‘‘Gaussian” here refers to the pulse shape as a function of time and not to 
the statistics of the pulses. 

The third type chosen is that of a square pulse envelope of finite duration. This 
exemplifies noise whose values are independent when separated by a sufficiently long, 
but finite time. The model fits the sonar and radar clutter problem when relative motion 
between the transceiver and scatterers is slight. 

Notice that, in all these cases, the semi-invariants can be expressed as functions, in 
fact powers, of the input covariance functiont. Accordingly, the time enters the output 
covariance only implicitly through the input covariance function. 

3.3 Derivation of the output correlation function. The total covariance function 
of the output is the ensemble average of the product of the outputs as two times separated 


*The alterations in the form of the covariance function which are necessary to take account of phase 





variations in time can be found in [11]. 

**H. Wallman has shown that the pass-band of cascaded networks whose individual step function 
responses have no overshoot tends toward the gaussian [17]. Furthermore, it seems likely that one can 
extend these results to any cascaded network except those tuned for a Butterworth response. 

+The linear model is one extreme of possible time behavior in that its semi-invariants of all orders 
decrease no faster than the correlation function. One may conjecture that no Az2(t), see Table 1, can 
decrease faster than the square of the covariance function, in which case the two other models represent 
the opposite extreme; however, we have been unable to prove this. 
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by an interval ¢. The average can be expressed as a suitable integral over the charac- 
teristic function. We have [14, 18] 


. : a = 
Rett) = @(Vi)g(V.)) = = I dé, dé, f(t) f(vé.)(exp (@ViE, + tV ok.) (24) 
tar ; 


ve 


sf 


where the average defines the characteristic function and 

fi) = | g(V)e*"* aV, Imé < 0 (25) 
in which g(V) is zero for V less than some V, , and of no greater than exponential order 
at infinity, and C is a straight line in the complex é-plane parallel to the real axis and 
lying below the singularities of f(7é). 

Since the input signal ensemble is the sum (see (21)) of sine wave and noise ensembles, 
which are independent of each other, the characteristic function of Eq. (24) becomes 
the product of the separate characteristic functions of the signal and of the noise. The 
results of Sec. II enable us to obtain the covariance function for non-gaussian noise, 
if we can obtain the corresponding one for non-stationary gaussian noise. Since the 
non-stationary gaussian case is only a very slight extension of previous work, it can 
easily be solved by methods briefly described below. 

The characteristic function of the sine wave signal Ay cos (wot + ¢) is [14] 


Fi(E, , & 5 Ds = Jo(Aoléi + & + 2k, cos wot]'””) (26) 
and the characteristic function of the noise is given by the exponential of Eq. (17). 


In order to find the zonal structure of the output, the two characteristic functions can 


be expanded in Fourier series, giving for the characteristic function of Eq. (24) 
7 V1 a + Wot 
Falb, bo 5) = exp( — Wert veke) yy ye OO eal Aoki a Ack) 
W7oe:E2) COS (Mwot + na) + Lj mn) (Wrobié2) COS (Mat — na)], (27) 
where the Neumann factor e,, equals one for m = 0 and is two for all other values. 


For the output of the receiver, only the low-frequency zone, m = 0, is required*. 
If now the modified Bessel function is expanded in a power series, the double integral 
of Eq. (24) separates into two single integrals, which have been evaluated previously [19], 
so that the covariance function of the low-frequency zone is 


ae 2 € (y'rs tia" Vivo\"”, =e v\~ 
tal) = y & Fy (k + n)!n? ( ) @p.)"( - es 


Sig 
iFy(—vp/2 +n + k3n + 1; —p) Fi(—~P/2 +n + kj n + 1; —prz), (28) 
where p; is the signal-to-noise power ratio, Ao/2y; , ,/, is Kummer’s form of the con- 
fluent hypergeometric function, and (a), = I'(a + n)/T(a). 


3.4 Strong signals. When the signal is strong, i.e. p; is large, Eq. (28) can be 


simplified by using the asymptotic expansion of the confluent hypergeometric function, 


*Note that, if the ordinary asymptotic form of the characteristic function of the noise is used, all 
the harmonics must be retained until the cross-terms between the frequency-dependent semi-invariants 


have been taken into account. 
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which gives 


Br’(y +1) [A5\’ =< € , 
R(t foo — a ( ) ae - es 2 = 2 ~ Jm\"t2k . 
WY) = 4T"(v/2 + 1) \4 > dX KV + nyt (~?/2)e(—¥/2)nsulto/p)"" C08 ne 


‘oF )(—v/2 +n +k, —v/2 + k; 1/p,)-2Fo(—v/2 +n + kh, —vp/2 +k; 1/p.), (29) 


where the formal hypergeometric function ./, is defined as 


n 


oF ,(a, 6; x) = + ks (a),,(b),, = (30) 
The asymptotic expansions used here are valid only when p > | —v/2 + +k |. How- 
ever, it can be shown [20] that the remainder contributed by the various possible series, 
e.g. when n > k > p, k > p = n, etc., is of higher order than the remainder indicated 
by taking a finite part of Eq. (31), provided that p > | —v/2 + + k | for all terms 
which are included. 

The differentiation of the gaussian result to obtain the non-gaussian one is readily 
accomplished with the help of the formula 


d 


dr F(a, b; ex) = abe .F,(a + 1, 6 + 1l;ex). (31) 
he 


When the leading terms of the power series are collected together, we find finally that 
the high and low density cases give the same result, viz: 


ery +1) / w\f , ‘ a : 
R(t) = ( ¥) 41 + 9/2)? LEM 4 (0/2) v/2)°1 + 1?) 
+ 1) \ P Pp 


iT (y/2 P 2 

+ (y/2 — 1)°(1 + 7, cos 2a) + 4(v/2)(v/2 — 1)ro cos a] (32) 
(y/2)°. : : 

4. SLAP, [(v/2 — 1)?Aso + 2v/2(v/2 — 1)[Asi(é) + Ais(0] 
3p ¥ 


+ [2(v/2)? + (v/2 — 1)*]Ao2(t)] + Op “*y »}. 


The order of the remainder is different as the noise density becomes large or small. The 
terms in p”’ are of order p* and p” y _'; the next set of terms are of order p*, p* y"', 
and p *y °. For nearly gaussian noise, y is large so that Eq. (32) can be used if p alone 
is large. For low-density noise, however, yp must be large as well as p. 

In physical terms, yp is the ratio of signal power to the noise power in a single noise 
pulse. When the noise density is small, there is no noise at all for an appreciable part 
of the time; thus the requirement on yp means that the signal must be relatively strong 
while the noise pulse is present, which is here a more stringent requirement than that 
it be strong on the average. 

When the signal is strong, there is, as usual, a modulation suppression effect [14] 
in that the leading term is the one for signal alone, with noise-entering as a first-order 
correction. Here, however, there is an additional suppression effect on the statistics of 
the noise. As Eq. (32) shows, the non-gaussian nature of the noise is first discernible 
in the second-order terms. 

3.5 Weak signals. The results when the signal is weak are more complicated because 
the low-density and high-density cases can no longer be described by the same formula. 
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If the confluent hypergeometric functions of Eq. (20) are expressed in series, one sum 
can be obtained in closed form, leaving a triple power series in the p; , viz; 


' Bry + 1) (yiv2\” ie Ce il 
R(t) = z= ( " } >> “kl dint n)! 


~~ 4T°(vp/2 + 1) f =o kod U0 


l+n/2 J,2 27/2 2.2 
ee (ee) COS Na 40, 0, 0; y 7 \ (33) 
: 1W2 1¥2 





where 

Fla, b,c; 2} = (—v/2)esnsel—v/2):anz0 
Fi(—-~e/2+k+n+a,—7/2+l+n+b);n+1+0¢;2). (34) 

Applying the differential operator of Eq. (17) one finds for the covariance function for 

nearly gaussian noise 


B°T?(v + 1) (4) > — = €,(—1)**! cos na _ 











R(t) = 777 (p/2 + 1) pp 2D ey l\na+h! n+ D! 
rip { 510, ,O;r}+= rk ao}, (35) 
where 
G(t) = AgoF{2,0, 0; 73} + [Asi(t) + Ais(0) [2s (1, 0, 0; 70} 
+ “1 §{2, 1, 1; 73} | + A, «| Me F{0, 0, 0; r5} (36) 
r eit 51,1, 0570} + + rt 2) 7t% 20} | 


The hypergeometric functions reduce to polynomials when » is an even integer and to 
complete elliptic integrals when » is an odd integer. 

When the noise density is small, one must use the characteristic function Eq. (20), 
which involves much the same sort of derivative as does the high-density case just 
treated. In addition, however, there are the time intervals in which no noise pulse or 
one noise pulse occurs. If there is no noise present, the correlation function is that for 
signal alone, which was found as the first term of Eq. (32). 

If a single noise pulse alone is present, the exact rl of the distributions must be used. 
On the other hand, when the signal is present at the same time as a single noise pulse, 
the covariance function is less sensitive to the statistics of the noise, and it becomes 
possible to approximate the noise distribution in the usual way. For a single pulse the 
covariance function is most easily found by not carrying out the average until the end. 
Thus the characteristic-function of the noise can be written as 


Fi(E, , &2 5 1) = (Jol ligt + v22 + Eww. cos (wot — a)]'”*), (37) 


where the average is over the parameters of v(t, — 7), an individual noise pulse. 
This characteristic function can be expanded in a Fourier series by the usual Bessel 


function addition theorem, whereupon the double integral for the low frequency zone 
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splits into the product of single integrals each of which is of Weber’s type [21], so that 
the correlation function is 





BT» + 1) (¥)’ (viv2)_ (38) 


41/2 + 1) \2/ 7@?y 


If now the differential operator of Eq. (20) is applied to Eq. (33), we obtain finally 
for a weak signal and low-density noise 


ri) = 2Pe+D (4) { pre” 
’ 41*(v/2 + 1) \2/ \(r*(v/2 + 1) 


4 i (vis) 


vy’ 'T°(v/2 + 1) ’)’ 

















c 2 @ k+ln 
4 €,(—1)*"" "ro cos na k+ltn 
> k=0 l=0 n! (n + k)! k! (n + l)! Th (39) 


ly" oly,» — n — k — DS{0, 0, 0; 73} 


1 l—rtnt+k+ 
+ 3p ¥ ‘oy,v—-l—-n—-—k— IG 


—iv—k-—l—-n@w-—k-l—-—k-—-1}) 


I 


—psnt+kel 2 
"y gy,y—-n—-—k— ps(0, 0, 0;r31}, 
where G(t) is defined in Eq. (36) and 


«o k 
¢(y, 8) = - k*e~” : (40) 
k=1 , 
y is a polynomial in y and exp (—y) for positive integral 8, and an iterated integral of 
the exponential integral for negative integral 8; for other values, a closed form is un- 
available. 

3.6 The quadratic detector. The output correlation function after a half-wave 
quadratic detector is particularly simple. When v equals 2, the output voltage is the same, 
except for a scale factor, as that from a full-wave squaring device, because the noise and 
signal distributions are symmetric about v = 0. The covariance function of all zones is 
thus a fourth-degree moment of the input. The covariance function of the low-frequency 


zone is then 


> yy? 2 2A22(t) 2 
Rif) = a 1 + To + 3y + 2p(1 + To) + Pp ’ (41) 
where p = A3/2y, the input signal-to-noise ratio, and rg = ro(t). The equation is exact 


for all values of signal-to-noise ratio or noise density. Finding the output covariance 
function as a moment of the input voltage is by far the simplest approach; however, the 
method may be applied only when » is an even integer. 

Results for the three noise models listed above in Sec. 3.2 are given in Table 2 and 
illustrated for the exponential and linear models in Figs. 2, 3, and 4. The figures show 
that the square-law rectifier reduces the amount of correlation compared to that present 
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TABLE 2 
Covariance functions for the quadratic detector. 


1. Exponential model: 
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Fic. 2. Output covariance for the quadratic detector exponential input covariance). 


in the input for all values of p and y. As p increases, the curves for all values of noise 
density increase too. This is because the signal X noise intermodulation terms become 
progressively more important as compared to the noise X noise products and, since the 
s X n contribution is more correlated*, r,.,.(¢) 1s increased. For large p a non-gaussian 
*The expression “more correlated” is used throughout the paper to express the fact that one nor- 


malized output correlation function is greater than another for the same value of ¢ when both are derived 


from the same input correlation function. 


1958 


RECTIFICATION OF NON-GAUSSIAN NOISE 
1.0 


409 
T 





0.8-- 


OUTPUT COVARIANCE 
QUADRATIC DETECTOR 


GAUSSIAN INPUT COVARIANCE] 
\ pt? 
0.6 





| 1 
0.5 1.0 





2.0 
B+ 















08 


OUTPUT COVARIANCE 
QUADRATIC DETECTOR = 
LINEAR INPUT COVARIANCE 


0.2-- 





1 
0.2 0.4 


Bt 
Fic . 








1. Output covariance for the quadratic detector (linear input covariance). 
suppression effect enters, so that the difference between gaussian and non-gaussian 
noise is less than at low signal levels. 


The effect of the non-normal nature of the noise here is to reduce the output correla- 
tion compared to the corresponding gaussian noise for the exponential and gaussian 
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noise models and to increase it for the linear model. Because of the simplicity of the 
expression for the covariance function, it is easy to see why this is so. For gaussian input 
noise, the quadratic detector output has one distorted noise term (ro , for n X n 
modulation). For non-normal noise, an additional term enters, containing the semi- 
invariant A,.(¢). This term adds to the ‘‘scrambled’’ noise for the exponential and gaussian 
models, and to the undistorted noise for the linear model. Clearly, the output correlation 
is usually decreased by non-gaussian statistics since, for nearly all variations, A,.(t) 
will decrease faster than r,(¢). The linear model is the extreme case in which this is not so. 

3.7 The linear detector. Obtaining the correlation function for a half-wave linear 
detector involves all the considerations, if not quite the full algebraic complexity, which 
attend the general case. No elementary method of obtaining the output covariance 
function is possible for a half-wave linear detector. However, when the general vth-law 
results of preceding sections are specialized to vy = 1, considerable simplification occurs. 
The hypergeometric functions of the general result reduce to complete elliptic integrals 
and the results for large and small density become much alike. Table 3 lists the cor- 
relation functions for the three noise models. The three cases of a strong signal, a weak 
signal with high-density noise, and a weak signal with low-density noise are listed sepa- 
rately. The expressions in the table are the leading terms of series in the noise density 
and signal-to-noise ratio p. The nearly gaussian and low-density results hold over a 
considerable range of variation of y, because the numerical coefficients of higher terms 
are small. The strong signal formulas of the table are reasonably accurate down to 
values of p of about 2. The weak signal correlation function, on the other hand, will 
hold only in a restricted range, for p no more than about 1/2, because only the first 
two terms of the ascending p series have been used. 


TABLE 3 
Covariance functions for the linear detector“* 


a. Strong signal, p > 1 








1. Exponential model: ro = exp(—B | t 
i.  2'%ypf. . ltn, 1 =| a F 9 | 
R i) = r it + 2p + 8p l—~ff% + 44 (1 To + 27) 
4! offi. — =" + ari) | + 0¢ -4 -»} 
6p" = ali 
Gaussian model m™=e- c 
Ri) = 260 +E, LU . — To)” + poll — 4r3” + 3rd | 

7 2p 8p" 


1 ~ = 3/2 m2 = -2 
+ — 6p [a — r5)(1 — Pr) +a —— 1 (5 — 8ro — 4r0° + vi) + 0(p", 7) 


*Here K(ro) and E(ro) are complete elliptic integrals of the first and second kind, respectively, and 
Biro) = E + (1 — ri)K, riC(ro) = (2 — ri)K — 2E, riD(ro) = K — E. (See Jahnke and Emde, 
“Tables of Functions,’’ Dover, New York, 1945, pp. 73-80.) Also ¢(y, — 1) is obtained from Eq. (40). 
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b. Weak signal, low-density noise, p < 1,7 <1 


1. Exponential model: ro = exp (—B | t}) 
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c. Weak signal, nearly gaussian noise, p < 1, y > 1 


] 


1. Exponential model: r, = exp (—B | t}) 


k(t) = P v ,E(r + rjB(ro) — (- oa [2K (r.) — K(r,)] 
T 3) 


37 
. ¥ Gaussian model: r — exp(— ES 
Bry l , = 
Ri ¥ dF ro) + r2B(r,) + - 5s [4r,D(ro.) — (A — ro) K (ro) — Aro D(ro)] 
T Sy(ar 
. . l , 5 ira. oe 
+ p| E(r B(ry) a {11 — Sry + Ary) E(r%) 
8y(1 i] 1) 
+ 7,(3 $7 + 7, )BO | +TUP,y ) 
) 
1— Bit até} <I 
3. Linear model: : 
0 Piet 4 
Ri) = 2¥ i e@,) + BG) — <2 [Ke.) — 27D )] 
Tv oY . 
E ; Bro) + . aati 7 3 Br D(r i + 0( a > 
Pp “7 T ToeD\lo) 7 84 es 7 0 P;Y7 


The normalized output correlation function is shown in Figs. 5, 6, and 7 with r...(0) 
for gaussian noise of the same input power shown in dashed lines. One may note, as 
for the quadratic detector, that increases in p cause the amount of correlation in the 
output to increase and cause the difference between gaussian and non-gaussian noise 
to decrease; however, the results have several features which are distinct from those 
for the quadratic detector. For example, as the noise density decreases, the leading 
term in the low-density series, (v’(t))v’(t) + ¢)), which for the linear detector is r,(t), 
becomes increasingly important, so that the output has more correlation than the output 
for a gaussian input noise for all noise models. This effect will occur as the noise density 
becomes sufficiently small for fixed signal power, even though the (average) signal-to- 
noise ratio is large. When yp becomes small, the weak signal series must be used. 

Determination of the extent to which the nearly normal and low-density results 
together cover the range of y is an important question which the covariance functions 
for the linear detector partially answer. Apparently the shape of the noise pulses has 
considerable effect. For gaussian pulses, the two forms of the correlation function are 
in close agreement (too close to be resolved in the figure); for the exponential and linear 
pulses, the agreement is not so complete. The figures tend to show that when the noise 
density is about 1, the two forms of the covariance function will be only roughly equal, 
but that no gross differences occur. 

3.8 Effects of non-gaussian statistics in general. The linear and quadratic detectors 
are the most important, so that the results applicable to them have been presented in 
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Fic. 5, Output covariance for the linear detector (exponential input covariance). 
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Fig. 7. Output covariance for the linear detector (linear input covariance). 


considerable detail. A qualitative look at the output correlation function for arbitrary 
v will indicate the over-all features of these results and indicate the place of the linear 


and quadratic detectors in the general scheme. 

Comparison of the covariance functions of gaussian and non-gaussian noise can be 
carried out for the weak signal cases by comparing the leading terms of the expansions 
of the hypergeometric functions. Separate expansions must be used when r; is near zero 
and when r, is near one. Since in both limiting regions, the same results emerge from 
the analysis, we may expect with a reasonable degree of confidence that the results 


t), and hence all ¢. If the signal is strong, all the functions involved 


hold for all values of r 
are elementary, so that the comparison is much easier to make. In both cases, the algebra 
lengthy, and not especially illuminating, so that it will 


is at once straightforward, 
22]. Figure 8 shows the results of such a comparison. In 


not be reproduced here [11, 
the shaded region the non-gaussian and gaussian correlation functions intersect each 
other, so that neither can be said to be larger, although for most values of ¢ the non- 
gaussian output is less correlated. No scale has been indicated on the vertical axis because 
the position of the boundaries depends on the density of the noise. Deviations from the 
gaussian result are small for large signals and large for small signals. 

3.2, the noise models chosen here represent nearly the extremes 


As mentioned in Sec 
of variation, so that we may expect the same general type of behavior for other shapes 
of noise pulses. Unless the semi-invariants of high degree have the same time structure 


[e.g. in the linear model, they all equal r,(¢)], the type of behavior shown by the gaussian 


or exponential noise model should be representative. 
3.9 The output spectrum. The variation of the output power as a function of 
frequency is an important statistic of the output. It can be obtained as the cosine trans- 
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Fig. 8. General effect of non-gaussian statistics. 


form of the covariance function, 


Wif) =4 [ " RY) conet at 


R(t) = [ W(f) cos wt df. (42) 


The covariance function is the quantity found directly from the analysis of the detector 
and therefore has been presented in greater detail than the spectrum will be. In general, 
the only tractable expressions for the covariance function are approximate, and the 
transform even of these can be taken only approximately. The covariance function 
and spectrum are equivalent in principle; and the propagation of error in taking cosine 
transforms can be minimized by giving results primarily in terms of the covariance 
function. 

The general behavior of the spectrum can be inferred from the reciprocal spreading 
property of Fourier transforms, i.e. the more a function is concentrated in one domain, 
the more it is dispersed in the transform domain [23]. 

A more quantitative measure of the spectrum can be found by the use of equivalent 
rectangular bandwidths. The low-frequency output spectrum has a maximum at zero 
frequency and decreases smoothly as the frequency increases. One can construct a 
spectrum having a constant value over a finite bandwidth equal to the zero frequency 
value of the actual spectrum and the value zero elsewhere, such that the total power 
in the actual and rectangular spectra is the same for both. Then the total bandwidth 
of the rectangular spectrum is the equivalent rectangular bandwidth of the actual 
spectrum. It is customary not to include the 6-function at zero frequency representing 
d-c power in the actual spectrum. 

Explicitly, one has 


[ wan ag 


fo _ ~ W,(0) ’ (43) 
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where W,(f) is the spectrum minus the d-c contribution. In terms of the correlation 
function, the bandwidth is 
R(O) — R(o~) 


} | [F(t) — R( 0 )] dt 


J0 


E [ ral at J (44) 


Table 4 gives the ratio of the output bandwidth to the input bandwidth for all the 
curves of Figs. 2 to 7. The entries in the table give a quantitative measure of the difference 


TABLE 4 


Spe ctral bandwidths*. 


Quadratic Detector 


Exponential Gaussian Linear 
covariance function covariance function covariance function 
p = 0, all y 2 00 p = 0, all y 1.41 p = 0, y= 1.50 
p= 5 = Ll dea p= a Y= me ey p= 0, Y= l 1.20 
p = .d, y= 1 1.50 p = .9o, y= 1 1.28 p =1, 7, = L.i2 
p=! y= 1.20 p=1, y = « eh p =1, y= 1 1.09 
p=1 y= 1 1.33 p=! = | 1.19 p = 2, y= « 1.07 
p=d y = 1.05 p=9d = 1.03 p = 2, y= 1 1.06 
p 5 y = 1 1.09 p = 5, y= 1 1.05 p = 5, y= o@ 1.03+ 
p=5 y= 1 1.03 — 
Linear Detector 
Exponential Gaussian Linear 
covariance function covariance function covariance function 
p=0 ‘ 2.25 p = 0, y = 1.48 p = 0, y,=o0 1.59 
p= 0 = 2NG 1.86 p= 0, y¥=2 1.42 p= 0, y,= ILD 1.28 
p = 0, y¥ = LD 1.55 p = 0, y¥ = 1.5 1.19 p=0, y=I1NG 1.18 
p=0 > l 1.25 p = YU, y= 1 1.15 p = 0, y= .0 1.08 
p = 0, y = .D 1.1] p = 0, y = .O 1.09 p= 0, y = 2 1.05 
p= 5 y= 1.44 p=.0 Y= 1.12 p=1, 7, = @ t.13 
p=. y¥=1 1.37 p=.9, Y= 1.08 p=1, y¥=I1 1.14 
p=2 Y= o Br. p = 2, y= 1.06 p = 2, y= @ 1.05 
p=2 7 ] 1.10 p=2, y=! 1.08 p=2, y=1 1.08 


*The entries in the table are the ratio of the bandwidth at the detector output to the bandwidth 


at the detector input. 
NG means “nearly gaussian”’; LD means “low density” 


between gaussian and non-gaussian noise. Since the bandwidth is larger as the covariance 

function is smaller, the results are just opposite to those of the previous section. 
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Appendix A 


Nearly normal distributions in terms of parametric derivatives. An expansion 
of a nearly normal distribution in a series of derivatives with respect to the second 
moments is useful when the noise ensemble has a phase uniformly distributed, inde- 
pendent of the envelope. As has been pointed out in Sec. 2.3, nearly normal narrow- 
band noise waves have such a structure. 

The terms of fixed degree, M = m + n, in the moment expansion of the second- 


order characteristic function (Eq. (3)) are 


_ R’ “ - n \/ se \m/s n 
S itt, cos” (wot) — 8;) cos” (wo[to + t] — 82))(ae,)"(7E2) 


“Mo (A.1) 
= ([2€,R, Cos (Wolo — 81) + 2&2 lk. COS WWo[to + t] — 6,)]"). 
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The average over wf) — 6; can be carried out as a contour integral around the unit 
= exp {7(wol) — 6,)}. There is one singularity inside 


circle, after the usual substitution, z 
the contour, a pole of order M + 1 at the origin; the integral is zero if M is odd, and, 
when / is even, gives 


M 
cos. @),,,. 2 ee os : ul 
Wn (((7t,R,)° + (t&.R>) + 2(2ER,)(t€.R2) cos wot — 6)]"’”). (A.2) 


For a non-random phase change, to which we here restrict ourselves, the expansion 


is found compactly through the formal use of tensor notation including the summation 


convention. 


Let the vector {X,}, (¢ = 1, 2, 3) be defined as 


(A.3a) 


iii . Re, . HR}, 
and let 
(mt = {£? | 2t.£, cos (wot — a), £3}, (A.3b) 
fa @) 4 
{D;} =(—-, —, ol (A.3¢) 
lay, Oro OY2) 
so that the expression (A.2) becomes, with M = 2p, 
‘—1)’, er (—1) 2p 7 re a 
-(cos’” 6(X,2;)") = D\0 cos” @)(X, a es 8 2 (A.4) 
(2.)! 


2p)! 
The first factors of (A.4) can be recognized as an envelope factor for the moments [from 
Eq. (12)] so that we can write 


(—1) = soit sats » 
———" in, °° Cons eM (A.5 
(2p)! 1 rs 
where 
Sin, @2 ae (cos” 6 6 a tie he ea a 
lisa 


Thus the moment expansion of the characteristic function, which in genera 


double series [see Eq. (3)], can be written as a single series when the noise is narrow band, 


: (—l)" a se , 
et) =e ee ee ae (A.6) 


3y comparing terms with the semi-invariant expansion of the characteristic function, 


one finds that, similarly, 


‘ — (—1)* = a e 
Fig, ,&;t) = expt > Qor alee ee a oe (A.7) 
where 
L,, ge ls Asq-zn; rni-g@, 


and, in particular, L;; = Ag-i-i.i+i-2 » Dige = Mo-i-g-e.isiee-s » (t, j,k = 1, 2, 3). 
If the exponential in Eq. (A.7) is expanded in its power series except for its first term, 
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it is readily seen that each =; may be replaced by a 2D, , giving the desired expansion 


whose initial terms are 
: J 9? 98 
Fé, ,&:3 = \! + A L;,;D;D; + 6! L,;,D,D,;D, 


9! | AS 
+} 3 LiL D.D,D.Di + oS exp {— 41,2.) | _ 


lIveeded 
‘Fe?r "re 


from which the manner of formation of further terms is evident. 
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SOLUTIONS OF THE HELMHOLTZ EQUATION FOR A CLASS OF 
NON-SEPARABLE CYLINDRICAL AND ROTATIONAL COORDINATE 
SYSTEMS* 

BY 
VAUGHAN H. WESTON 


University of Toronto 


Abstract. In cylindrical and rotational coordinate systems, one of the variables 
can be separated out of the Helmholtz equation, leaving a second order partial differential 
equation in two variables. For a class of the coordinate systems, this equation is reducible 
to a recurrence set of ordinary differential equations in one variable, which are solvable 
by ordinary methods 

1. Introduction. The usual method of solving the Helmholtz equation 


Vyvy+ky=0 (1) 


separated into three ordinary differential equations, each of which can be solved. How- 
ever, it is shown by Eisenhart [1] that the method of separation of variables in euclidean 
space is applicable to only eleven coordinate systems, generated by confocal quadrics 
or their degenerate forms. Exact solutions of the Helmholtz equation have thus hitherto 


in three dimensions is the method of separation of variables, in which the equation is 


been limited to separable coordinate systems. 

In this article a method is presented for solving (1) for a class of non-separable 
cylindrical and rotational coordinate systems. 

2. Cylindrical and rotational coordinate systems. The definitions of cylindrical 
and rotational coordinate systems are those given by P. Moon and D. Spencer [2]. 
Let z = F(w), where z x, + ty, , and F is any analytic function. Separating real and 
imaginary parts, one obtains 


£,(u, v). (2) 


The curves u = constant, and v = constant give rise to an orthogonal family of curves 
in the z-plane. For each function F, a cylindrical coordinate system and one or two 


oordinate systems can be formed. 


rotation ] 
The cylindrical system (uu; , %2 , Uz) is given by the equations 


| é l; u 
y= £(uU Us (3) 
g Us y 
and the rotational coordinate systems (u; , U2, Us = ) by 
1 E\(U, , U2) COS 9, 
y = &,(u, , U2) sind, (4) 
2 E(u, Us) 








*Received December 20, 1956; revised manuscript received April 8, 1967. 
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and 


I 


z = £,(u, , U2) Cos, 
y = £,(u, , U2) sing, (5) 
E,(u; , U2). 


System (4) is found by rotation about the y,-axis and system (5) about the 2,-axis. 
The z-axis is then taken as the axis of rotation. The metric coefficients for the cylindrical 
coordinate systems have the property that 


_— 2) (28) _ , 
h; = hs — (2 = Dp ’ hz = l (6) 


and for the rotational coordinate system (4) 


i= (2: + OU» 


Zz 


(7) 
hs = &:(uy , Up) 
and rotational coordinate system (5) 
: : dk, \” dg, \* 
Ou, Ou 
(8) 
hs = &2(uy , U2). 
Also since F is an analytic function £;(u, , v2) has the property 
i, 4 Ho, 
du Ouz 
, . (9) 
0& _ Oks O  —__—OEa 
Ou, OU, ’ OU Ou, 
Now rotational and cylindrical coordinates have the property that one of the variables 


@ and uz respectively) can be separated out of (1) leaving a second order partial differ- 
ential equation in the two variables u, and u, . We shall show that there exists a class of 
coordinate systems for which this equation can be reduced to a recurrence set of ordinary 
differential equations in one variable. 

3. Solutions of the Helmholtz equation for rotational coordinate systems. 
Theorem I. Tf for a rotational coordinate system there is a u,;(¢ # 3) such that 

h;(u; , U2) has the property that 


Ih M, 5 : , 
“<2 = DY f,(u,){halur , us)]’, (10) 
Ou; seh; 
where NV, , M, are integers, and N, < M, and if the metric coefficient hj has the form 
Ms 
hi = > c,(us){ha(u , u2)]’, (11) 
s=No 


where N, , M, are integers and N, < M, then 
Wu, , U2, d) =e > a,(u,)[hs(u; , u2)]" (12) 


r 
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is a solution of the equation V*y + k*y = 0, where a,(u;) satisfies the set of recurrence 


relations 


da r+1-—M, da r+2—-M, 
a i ° 
© ——! (2g + I)f-+1-e(us) + =.) 
au; q=r+1-N; du; q=r+2-N, ‘ 
(13) 
r—M, 
Ag(UdCr+2-Us) +k? DS a,(uie,-(u,) = 0. 
a=r-N; 


Proof. The operator Y° in rotational coordinate systems has the form 


I 0 oO 0 es hy x 7 
hihs 152 (h. =) + OU2 (h. 2.) +r hs a} 


3 


Hence using the expression for y given by (12), Eq. (1) reduces to 


( ‘ , 2 
| ] : (n, c += (i, 2) — te + ih. > a,(u)(hs)’ = 0. (14) 
\ 3 / 


h7h; | Ou; * OU, OUz OUg/ 
Multiplying (14) by hj exp(—zu@) one obtains the following 


fF af a2 o a e ¢ 
Jo a , 1 dhs 0 1 oh; O oh 


- 4 - -— 1 4 en? -a,(u,)\(hs)’ = 0, 
> (Ou; OUz = hg OU OU hz OUs2 OUz # i if (usd) ; 


d’ la, ; ; .* a ie 
¥ ay SE + 2 ny Se + (in: — 4 Ms a)” 
( du; du; Ou; hg, 


yi 9h. \” ‘ah; \" yh, hs 
+ rhs) ‘a ( = (2) += »(%) + hy 3? + hs * 4 
\OU; OU2 Ou, U2 


h,)" a, | dhz dhs ; dh; dhs ' (h3)" dh; da, 0 
+ T\Nsz) = a a ae ls “¥ a = . 
hz Lu, Ou, OUy OUs h, ou; du; 


Using relations (7), (8) and (9) one obtains 


( > 
,da, da, --1 Ohs ein 
> §(hs)” at + (2r + 1) = (hs) * = + r’a,(hs) “hi 


du; au, ~~ Ou - (15) 
+ az —p Matt)" = Q. 


Use relations (10) and (11), and arrange (15) in power series in h3(u; , U2) to obtain 


vu 94+9— 


: s (h2)"S da, + ba da, (2g + I)fe+1-uis) + >. 


bye (q° — p)a,(U;)C,42—-¢(Ui) 
du; ss K du; Y KM )U, , 2-—¢ 


f, 
r N2 


\ 


r-M, 


+ pe k*a,(u,)c,-,(u;) ? = 0. 
q=r-N; 


Hence equate coefficients of (h3)" to zero. A recurrence set of ordinary, differential equa- 


tions is then obtained relating the functions a,(u;) 


2 r+1-Mi r—-M, 
da, + 2. _ (2q + I)fes: (u;) +k > ® a,(U,)C,—-g(Ui) 
du; Pp du; é a=r-N;z 
2- Ms (16) 


+ DY (¢ — wa, (uie,+2-(ui) = 9. 
q=r+2-N3 
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Thus provided that the ordinary differential equations given by (16) can be solved for 
the a,(u,;), the Helmholtz equation has the solution 


W(ur , U2 ,%) =e Do aus) [haus , u2)]”. (17) 


Here the summation over r represents a power series in the function (h;). 

There are two cases of practical importance, (i) lower termination of series (17) when 

M, > N, > 2and M, > N, > 1 and (ii) upper termination when N, < M, < 0 and 
NS ay = 4, 
For the case of lower termination the differential equation (16) is an inhomogeneous 
equation, the homogeneous part involving a, and the inhomogeneous involving 
d,-,; , @,-2 , *** ie. terms a, such that n < r. Thus, if a,_; , a,-., --+ are known then 
a, can be found. Now there is some number p such that a,_, = a,-2 = a,-3; = 0. Hence 
the inhomogeneous portion of (16) for r = p, is zero. Thus a, is a solution of the homo- 
geneous equation and can be found. Since Eq. (16) for r = p + 1 involves only a, and 
a,,, and a, is now known, then a,,, can be found. Hence if a, , a,,,; , «- + , a,-, are known, 
then a, can be found. 

Hence one has a solution of the form 


1 


ae 7 
Vi(u, , U2 ,o) =e” ys a,(u;)(h3) | 
r=p,p+1 (18) 
= e***(hs)” >» ay.n(hs)”. | 
N=0 } 
Similarly it can . shown for the case of upper termination that if g is such that a,,, = 
Berg @ Boas , the expression (17) becomes 
r=a,a-1 } 
ilu ur ,g) =e D> a,(us)(hs)’ | 
7 ] (19) 
| 


e'”*(hs) ‘= A_nv+(t;)(hs)~ 


The numbers p and g are determined from boundary conditions. In addition to the 
solutions of the form (17) there may exist other solutions described in the following 
theorem. 

Theorem II. If for a rotational coordinate system there is a u,(¢ # 3) such that 
h3(u; , Ue) and h? have the properties given by (10) and (11) respectively, and if there 
is a function B(u;) where t ¥ 7 ¥ 3; such that 


dB oh; _ — 8 9 
de: Oe = Biu;) > d,(u,)(hs)’, (20) 


where NV, and M, are integers and N,; < M; and 


d’ ‘B — 
= Blu;) > e,(u,)(hs)’ (21) 
ae P=. i 


where NV, and M, are integers and N, < M, then 
Wu , U2, ) = e'™*B(u;) Do b,(us)[hs(ur , us)]" (22) 
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is a solution of the Helmholtz equation V°y -+ ky = O where b,(u;) must satisfy the 


+ 


set of recurrence relations 


db, =~” @o ’ M 
— a > Fa + Dfssi-et 2 ( — w)Derr2-< 


du- ft du 
ul le r—Ms+1 
] p b.c, + >. Do oT ® b,(2g a 1) d ; Q). 


The prool is similar to that of Theorem 1. As before, there are two cases of practical 
importance for (22); upper termination and lower termination of the series. 
When M, > N aM, > NN. > 1, MM, > Ne > Cand N, > N,.=> 1 there.are 


solutions of the Helmholtz equation of the form 


and when N, < M, < 0,N, <M, < 1,N,< My, < O and N; < WM; < 1 there are 
solutions of the form 


a ee eo ek Sn) 24) 


Ti 


) 


The second type of solutions given by (23) and (24) are essential if h; is an even function 
of the variable u; . Then the first type of solution given by (17) are even functions of 4; , 
and if odd solutions in the variable u; are required, then this B(uw;) is chosen such that 
B(u;) is an odd function of wu; . 

4. Solutions of the Helmholtz equation for cylindrical coordinate systems. In a 
similar manner as fol the rotational case, solutions of the Helmholtz equation ean be 
obtained and are stated in the following theorem (where to simplify presentation the 
two cases corresponding to Theorem I and Theorem IT are combined): 

Theorem III. If for a cylindrical coordinate system there is a u;(¢ # 3) and &,(u, , v2) 


(k = 3) such that 


0& 
== = > b,(u lEz(U, y U2 | : 25) 


Ou a 
: va 


where N, and MW, are integers and N, < M, and if the metric coefficient h, has the prop- 


erty that 


where V, and M, are integers and N, < M, and there exists a function B(u;) where 


1+ j > 3 such that 


d-B We) = 
B 5 = Blu ) S é\u )[E.(uy » Uo) | 27) 
du dnd silos 
and 
dB Of ~~ 
— —— am See pd d,(u,) [E(u , U2) ]", 28) 
du, ou 
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where N; , AJ; , N, and M, are integers and N, < M;,;,N, < M,. Then 
V(uy , U2 , Us) = e™™Blu;) >> a(ud leu , v2)]" (29) 


r 


is a solution of the Helmholtz equation V*y + k*y = 0 and a,(u;) satisfies the set of 


recurrence relations 


da ia = da. , ree ; 
+ > r (2q)bu-u) + DO gq — Da,(uie,+2-(u,) 
au, Q 1-N au; @=r+2-N; 
r-M;, r—-Ms 
+(P —w) DY aure,-(u) + >> a(uve,_(us) (30) 
q=r-N, q=r-Ne 
r+1-M, 
—- z 2qa,(u;) d,-g41 = 0. 
qa 1-Ns 


A special case of the expression given by (29) corresponding to (17) for rotational coordi- 
nates, is given when B(u;) = 1 and d,(u;) = 0 fors = N, , N; + 1, --- , M; and 
e,(u;) = Ofors = N,,N,+1,---,I,. There are two cases of practical importance: 
(i) Lower termination, 7, > N, > 2,M, >N, >1,M,>N,>0,M3; >N;32> 1 
(ii) Upper termination, NV, < M, <0,N, << M,<1,N,S5M,<0,N; © M; < 1. 


For lower termination the expression given by (29) reduces to 
Wi(u; , U2, Us) = e'*Blu)) >> a(ud [E(u , u2))’ (31) 


and for upper termination 


Wi(u, , U2 , U3) = e'””*Biu,) > a,(u;)[E(ur , U2) ]. (32) 


5. Particular coordinate systems. The question of practicability of the above 
method of solving Eq. (1) for the particular non-separable coordinate systems under 
consideration is answered by the application of the method to the toroidal coordinate 
system. The author has obtained a complete set of solutions of Eq. (1) in toroidal coor- 
dinates, which satisfy the radiation condition and possess a ring singularity [3]. Besides 
toroidal coordinates, Eq. (1) can be solved for other well-known coordinate systems, as 
is shown in Table 1. The solutions given there, are independent of prescribed boundary 
conditions. The expressions given are in terms of power series of f(u;) hs(w , U2) and 
g(u;) &(u, , U2) instead of h3(u, , v2) and &,(u, , u,) according to whether the system is 
rotational or cylindrical. An appropriate choice of f(u;) or g(u,) simplifies the differential 
equations (13) and (30), transforming the homogeneous part into a recognizable form. 

Acknowledgement. The author wishes to thank the National Research Council 
of Canada for financial assistance through scholarships. 
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—NOTES— 


ON LINEAR PERTURBATIONS* 
By AUREL WINTNER (The Johns Hopkins University) 


If f(t) and g(é) are continuous functions for large positive ¢, and if the solutions 


of the differential equation 
a” + tl 2 =a ( (1) 


are known, how “small” must be the difference f(t) — g(t) (for large ¢) in order that the 
general solution of the differential equation 


/ 


y’’ + gidy = 0 (2) 
can be guaranteed to have the same asymptotic behavior as the general solution of 
(1)? The literature consulted answers this question only under assumptions which 
restrict (1) by certain conditions of stability.‘ No such assumptions are made in the 
following theorem (which, under stability assumptions, reduces to known results 


With reference to the coefficient function f(t) and a pair x = u(t), x = v(t) of linearly 
indepe ndent solutions of (1). let the coe fhicie ni function g(t) of (2) satisfy the following 


condition ; 


| f-glQulP+lof)dt < . (3) 
Then every solution 7 j(t) of (2) zs of the form 
y(t c,u(t) + c.v(t) + o(| u(t) | + | o(2) |), 1) 


where ¢c,; , C2 are integration constants (which can be chosen arbitrarily). In addition, the 


asymptotic relation (4) remains true on differentiation, that is 
y(t) = cu’ (t) + ¢v’(t) + o(| u’(t) | + | v’(t) |), )) 


where ' = d/dt. The o symbol in A(t) = j(t) + o(| k(t means that k ¥ 0 for large ¢ 
and that (h — 7)/k ~Oasit— o. 

The relation (4) reduces toy ~ c,u + c.vast— @ if (1) is stable in the sense 
z(t)| < ©, where ¢t > o, holds for all solutions of (1). In this case, condition 


/ |f—g|ldt<o@. (6) 


But the converse conclusion cannot be made (not even lim sup | x(t)| < © is assumed 


that 


lim sup 
3) is certainly satisfied if 


for z = uand x = v), since, when u(t) and v(t) happen to be “small” (for large ¢), then 
condition (3) requires of the “size’”’ of the perturbation f — g substantially less than 


what is required by the standard assumption (6). 
*Received October 11, 1956. 
1In this regard, see A. Wintner, Quart. App!. Math. 13, 192-195 (Secs. 2 and 5) (1955). 
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What is more, the general theorem applies without any stability restriction of the 
given problem (1). In fact, no matter what the coefficient and the general solution of 
(1) (that is, the function f and two, linearly independent, solutions u, v) may be, the 
perturbation f — g of (2) on (1) can always be chosen so small as to satisfy condition 
3), that is, so as to bring (2) within the range of the theorem. 

The proof of the theorem will consist of two steps. 

First, recourse will be had to an elementary lemma which goes back to Bécher* and 
which, in the binary case at hand, runs as follows: If the coefficients of a homogeneous, 


linear differential system 


p’ = a;(t)p + a,,(d)q, q’ = a2,(t)p + an2(t)q (7) 
are, for large positive ¢, continuous functions satisfying 
| |a,(t)|dt< o, where t=1,2 and k = 1,2, (8) 


then, corresponding to any pair c, , ¢. of integration constants, the system (7) possesses 


a unique solution (p, q) satisfying 
pthc, , gj—-c, as to @, (9) 


Next, the following rule of Lagrangian “variation of constants’ ( a rule which, 
being purely formal in nature, requires only the continuity of f(¢) and g(¢) on a /-interval) 
will be needed.* Let two solutions, x = u and x = », of (1) be so chosen that 
their Wronskian u(t)v’(t) — v(t)u’(t) (which is always a non-vanishing constant) becomes 
the constant 1, and define, in terms of the difference of the coefficient functions of (1) 
and (2), a binary matrix function |! a,,(¢) || as follows: 


a, a2) _ iy _ of mw —o\ (10) 


2 
Gs es u Uv 


Then y(t) is a solution of (2) if and only if there belongs to it a solution (p, q) of the 


case (10) of (7) in such a way that 
y(t) = u(t)p(t) + v(t)q(t) (11) 


becomes an identity. 
In addition, this transformation of (2) into (1) is a “contact transformation,” in the 
sense that the following differentiation rule holds for (11): 


y(t) = w'(Op + vo’ (Og (12) 


(in other words, wp’ + vq’ vanishes for all ¢)*. 
In order to combine this Lagrangian rule with Bécher’s lemma, note that the case 
(10) of the four conditions (8) is equivalent to the three conditions 


| |\f—g||w|dt<o, where w=u’,w,v’, 


For references, and for certain refinements, see A. Wintner, Am. J. Math. 76, 183-190 (1954). 

For a verification (and for a similar application) of this rule, see A. Wintner, Am. J. Math. 69, 
262-263 (1947). 

‘See formula (34) in the preceding reference’. 
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and that, since 2| w| S| u|* + |v |’, the latter three conditions are equivalent to the 
single condition (3). Accordingly, (3) assures the validity of the limit relations (9) for 
the general solution (p, q) of the case (10) of (7). But (11) reduces to (4), and (11) to 
(5), by virtue of (9). 

This completes the proof of the theorem. Its assumption (3) is independent of the 
choice of the two, linearly independent, solutions x = u(t), «x = v(t) of (1) which occur 
in (3). For, on the one hand, u(¢) and »(¢) cannot vanish at the same ¢ and, on the other 
hand, the ratio of | u*(¢#) |* + | v*(#) |? to| w(t) |? + | v(é) |? stays between two positive 
constant bounds ast > ©. This is clear from the fact that u*(t) = au(t) + bv(t) and v*(t) 
= cu(t) + dv(t), where a, b, c, d are constants of non-vanishing determinant ad — be. 

A particular, but interesting, case of the theorem results if condition (3) is strength- 
ened so as to make possible the elimination of the solution pair (u, v) occurring in (3). 
Such an elimination is made possible by using a well-known estimate, which was applied 
in more general forms by Liapounoff’ and others (L. Schlesinger, G. D. Birkhoff, O. 
Perron)” and which, when particularized to the case of (1), states that 


| a(t) | < const. exp tal | f(s) — 1 | ds} (13) 


holds for every solution 2z(¢) of (1). What then results avoids the implicit hypothesis 
of the theorem, namely, that (1) has already been solved. In fact, the resulting corollary 
of the theorem can be formulated as follows. 


If the coefficient function g(t) of (2) is so “close” (for large t) to the coefficient function 
f(t) of (1) that 


| g(t) — f(t) | exp ‘| | f(s) — 1| ds} dt < @, (14) 


then the asymptotic behavior of all solutions y(t) of (2) and of their derivatives y'(t) is given 
by (4) and (5), where c; , C. is an arbitrary pair of integration constants and u(t), v(t) is a 
pair of linearly independent solutions x(t) of (1). 


In fact, if (13) is applied to x = u and x = », then (3) reduces to (14). 


A LEBEDEV TRANSFORM AND THE “BAFFLE” PROBLEM* 
By C. P. WELLS AND A. LEITNER (Michigan State University) 


1. Introduction. This note is concerned with the application of the Lebedev trans- 
form to what we term the “baffle problem,” i.e. the problem of sound radiated by a 
vibrating circular disk in an infinite rigid baffle. The solution of this problem is not 
new. It has been solved by Sommerfeld [1] in terms of cylindrical waves, using the 
Hankel transform, and by Bouwkamp [2] and others in terms of spheroidal waves, using 
series representation. However, the Lebedev transform offers an equally straightforward 
method, representing the radiation in terms of spherical waves, and moreover is in 





5E. Picard, Traité d’ Analyse, 3rd ed., vol. 3, 1928, p. 385. 
*For the particular case (13) of the general theorem, see N. Levinson, Duke Math. J. 8, 2-3( 1941). 
*Received November 15, 1956. Sponsored by Office of Ordnance Research, U. S. Army. 
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itself of considerable inherent interest. It has been applied by Kontorovich and Lebedev 
[3] and by Oberhettinger [4] to problems of diffraction by a wedge, and by Leitner and 
Wells [5] to the problem of a freely vibrating disk. 

We use spherical coordinates r, 6, ¢ with the baffle in the plane 6 = 2/2 and seek 
solutions, independent of ¢, of 


V?ut+ ku =0, (1) 


where k = 27/\, \ = wave length, and u = u(r, 6), the velocity potential. The boundary 
conditions on u are: 


Ou _ ‘° a constant, whenr < a, 0 = /2, 
on 0, r>a,60=7/2 
where 0/dn is the normal derivative, together with the radiation condition r(iku + du/dr) 
— OQas7r— o~. We now try to represent u(r, 0) as: 


u(r, &) = i] u9(u)P_+, (cos 6) J,(kr) dp, (2) 


where P,,.,,(cos @) is the Legendre function, J,(kr) the Bessel function, and g(u) an 
unknown function of the complex variable 1» = o + ir. L is a contour in a strip of finite 
width surrounding the imaginary yp axis from ¢ — im too + i . The function g(u) 
is to be determined by applying the boundary conditions to (2) and using the theorem 
of Kontorovich and Lebedev [3] which states that if 


(lr) = ff wa(we'™*J,(07) du, (3a) 
L 
then 

riA(u) = [ o(krje**"’? HH (kr) dr/r. (3b) 


The conditions for validity of this theorem and other details can be found in the reference 
given. 

One now applies the boundary conditions to (2) only to find that the resulting A(u) 
is such that the conditions of the theorem are not satisfied and the integrals diverge 
However, there is an analogous theorem [6] in terms of modified Bessel functions of 
real argument which imposes considerably milder restrictions on A(u). This suggests 
making the transition k = —ty, y real and positive, constructing a solution using the 
modified functions J,(yr) and K,(yr) and then returning to real k and the original 
Bessel functions. Obviously the integrals would still diverge if the contour LZ remains 
unchanged, but they will converge if L is first deformed to surround the positive real 
u axis. This was recently demonstrated by Oberhettinger [4] and verified again by the 
result of the present problem. 

The representation of u(r, 6) corresponding to real y is 


Hur, ) = | wglw)P-4.» (cos 6) I,(07) dy. (4) 


The ‘baffle problem” can now be solved by means of the following Lebedev [6] 
theorem: If 
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g(r) = wA(u)I,(yr) du, (5a) 


then 


miA(u) = | g(yr)K, (yr) dr/r. (5b) 
The conditions for validity are that A(u) be an even function of yu, analytic in a strip 
of finite width including the imaginary axis and of decay at least as fast as 
r|**exp(— 2|7/1|/2), where e€>0O. 
2. Some properties of /,(yr) and K,(yr) as functions of ». The functions /,(y7 
and K,(yr) with real argument yr are entire functions with an infinite number of simple 


zeros [7]. For | u/yr | > |, 


Tyr) = [(yr/2)*/T0 + wl + O(u7")). (6) 
Hence J,(yr) decays rapidly as Re p» — + o. On the left half plane, as Re p > —@, 
it has T function-like growth since |/I'(1 + yz) — sin muI'(— pu)/z. Along the imaginary 
axis u = ir, I,(yr) has exponential growth as r > + @. 


The function K,(yr) is defined by 


K, (yr) = (2/2)[I_.yr) — L,Cyr))/sin rz. (7) 

As Re p — -&, I_,(yr) is dominant and as Re up > — o, J,(yr) is dominant. Hence 
K,,(yr) has T function-like growth on both right and left half planes. Along the imaginary 
u 2 “ 
axis K,(yr) decays exponentially as 7 — + o. Asymptotic forms for J,(yr) and hence 


for K,(yr) can, of course, be found from (6) using Stirling’s formula for (1 + xz). 
3. Solution of the ‘‘baffle” problem. We represent u(r,6) by means of (4) and enforce 
the bour dary conditions as modified by transition to real 6 a 


eo 3/2. ) . 
r<a, r’'v ; ‘ 
Pa a: 0 can 
Here A(u) = g(u)P’;.,(0), where the prime indicates differentiation with respect to the 
argument cos @. The formal solution is given by the inversion integral (5b 
awiA(u) = 1 | rK yr) dr. (9) 


The integral defines an even function of u and converges if yu is restricted to a strip of 
half-width 3/2 about the imaginary axis. Within this strip A(u) is analytic and 
.| A(z) | my eo '9'? 1797 a3 \ 7 — &. From this it is seen that A(u) satisfies the conditions 
of the Levedev theorem. 

The integral (9) can be expressed in terms of the Lommel functions and they in 
terms of their series representation [8]. Thus one can continue A(x) into the remainder 
of the u plane where it is analytic except at the points y = + (2n + 3/2), = 0, 1, 2, 

, where it has simple poles with residues [2”” v(—)"i/wy*”’] [[(n + **)/n!]. For 
Re » — + -&, we find also that A(u) ~ I.,(ya)/(3/2 + wu) sin mu. We now have A(u) 
defined for the entire » plane and are ready to convert the integral (4) into an eigen- 


function expansion. 
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Substituting for g(u) in (4) gives 


rue, ) = ff way SBS rn) a (10) 
in which pA(u) I,(yr) = (r/a)*/p for large uw on the right half Long Further the ratio 
of the Legendre functions behaves like e“’~*’”* when @ < 7/2 the only values of @ in 
our problem. Hence for r < a the contour I may be closed on the right and the integral 
is unchanged in value. The poles of the integrand are at 2n + 1/2 and 2n + 3/2 and the 
series of residues is the eigenfunction representation for the region r < a. 

For r > a the contour of (10) in its present form can not be closed. However the 


decomposition 
A(u) = (w/2)[A(— ») — A(u)]/sin ay 


allows the contour to be closed. Then 
miX(u) = v | rd,(yr) dr, (11) 
0 


and, by arguments similar to those used above for A(u) we find A(z) to have poles at 
—(2n + 3/2), analytic in the strip with growth like J,(ya)/(3/2 — wu), for large u. 
We now have 


Pale, 8) = funy) hee Kalan) du, (12) 


where pA(u) K,(yr) & uw (v/a)™, | m | — o, Hence the contour can be closed on the 
right for r > a. The poles of the integrand lie at 2n + 1/2, whose residues lead to the 
appropriate series of eigenfunctions in the space r > a. One sees now how the transition 
from integral representation to series of appropriate eigenfunctions resolves itself. When 
the transition from real y to real k is made we see that the expansion will be in terms 
of Bessel functions for r < a, in terms of Hankel functions for r > a. 

For completeness we record the eigenfunction expansion: 


r<a: ulr, 0=) = > a,Pons: (COS O)jonsi(kr) + >> b,Pon (COS 6)jon(kr), 
n=0 n=0 


C) 


r>a: u(r, 6) = +3 CrP on (cos 0) hs?’ (kr), 


n=0 


where j and h stand for the spherical Bessel and Hankel functions and 


a, = (2v/k)(2n + 3/2)(— 1)", 
b, = (2)*vai(2n + 4)(— 1)" [Tn + 4)/n!JW {He 4(ka), 83 ,2.44(ka)}, 


ka 


c, = ((2)tiv/k)(Qn + 2(— ITH + p/nty [ 2*Jons4(x) dz, 


0 
where W stands for Wronskian and 8} ,2,,; is a Lommel function [8]. Note that for r > a 
only the Legendre polynomials even in cos 6 appear, as required by the boundary con- 


dition. 
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GEOMETRIC INTERPRETATION FOR THE RECIPROCAL 
DEFORMATION TENSORS* 


By C. TRUESDELL (Indiana University and Universita di Bologna) 


In a finite deformation x = x(X), changes of infinitesimal lengths may be measured 
by the tensor C, where 
k m ry 7K >M k 
ds’ = Jim dx” dz™ = Crm dx” dX ; Cru = Dem@ 0 P (1) 


or by the dual tensor c satisfying formulae that follow by systematic interchange of 
majuscules and minuscules. Geometric interpretations of C and c have been given by 
Cauchy and others. In 1894 Finger introduced the reciprocal tensors C~* and c™', and 
recent exact work on isotropic elastic bodies employs them often. While formulae such as 


(CORN = g'"XAX (2) 


for their expression and use are known, geometric interpretation has been lacking. 
As is known, the correspondence between elements of area is given by da” = 2° x 2" 


dA*™”. where dA*™ is connected with the usual vector element of area dAx by dAx = 


‘ MP y 1/2 . 1 
exupdA » CxmpP (det Gor) €KMP:. Hence 
,7\2 : F PG RS 
(da) = Creal X” pt “ox' 2X’ s dA ? dA ; 
det q km KRS +r .3 UPQ 4 ‘ 
— 7 + one ( + €, ra€ XL rX.s)\ de npa€ , "zr px %G) dArx dA My 
ae Trt . 
(3) 


det g Ox’, x7, 2°) |? a 
~ det . ar, x. x gp XuXim dx dAn , 
A TUYV \h yy <b A 


Ke 


= det fe, gt ff _ . dAx dA xy . 
: ; ; ; - 
Comparing this result with (1) shows that the tensor C-'/det C™” measures changes 
of the magnitudes of infinitesimal areas in precisely the same way as C measures changes 


of infinitesimal lengths. 
A known principle of duality, which may be called the first principle of duality, 


*Received Dec. 4, 1956. This work was done under a National Science Foundation grant to Indiana 
University. 

1A formula which is essentially the next to last step in (3) was given by Tonolo, Rend. sem. mat. 
Padova 14, 43-117 (1943), Sec. V. 4, but he did not mention any connection with C™. 
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enables us to interchange the roles of X and x and thus obtain an analogous interpre- 


tation for c ‘/det c™' 


These results are equivalent to a second principle of duality: Any proposition on 
changes of length expressed in terms of C and ¢ yields a theorem on changes of area if 
C, c, and “length” be replaced by C~'/det C~’, c7*/det c™’, and “area’’, respectively. 

Of the many theorems that may be derived in this way, I record only one: The 
elements of area suffering extremal changes are normal to the principal directions of 
strain, and the greatest (least) change of area occurs in the plane normal to the axis 
of least (greatest) stretch; in fact, if the principal stretches dx/dX satisfy 4, 2 Az 2 As 
the corresponding ratios da/dA satisfy \.A43 S AsA,; S A,A2 . While this theorem is geo- 
metrically plausible, the first part does not seem obvious. 


BOOK REVIEWS 


(Continued from’p. 394) 


The theory of linear antennas. By Ronald W. P. King. Harvard University Press, Cam- 
bridge, 1956. xxi + 944 pp. $20.00. 

The basic material in this treatise on the linear antenna is founded on a graduate course given by 
the author, and also includes both the theoretical and experimental data of many other workers in the 
antenna field. Although much of the information has been published previously in various technical 
journals, the detailed and extensive compilation and evaluation of so much of this work may be regarded 
as a worthwhile contribution to the antenna specialist. It is the belief of this reviewer, however, that 
the general usefulness might have been greatly enhanced if the material had been separated into two or 
more books. 

The stated purpose of this book, which is basically mathematical in its approach, is to provide a 
bridge from the mathematician to the practical antenna engineer. An introduction summarizing the 
highlights in the historical development of the linear antenna theory is followed by a short chapter on 
the essentials of electromagnetic theory. 

The mathematical difficulties in treating the behavior of a single linear radiator as end-load for a 
two-wire line are considered in chapter II. Only antennas having a cross section which is small compared 
with the wavelength are investigated in this text. Particular attention is paid to the end effect compli- 
cations and cross-coupling between the line and the antenna. An approximate method of compensating 
for these effects is developed which employs appropriate lumped reactive elements at the junction 
between the line and the antenna. The characteristics of the isolated antenna are then studied in detail 


using several different formulations of the problem. The antenna impedance and admittance calcu- 
lations are presented in a number of different types of graphical plots as well as in useful tables of 
numerical values. Unfortunately, in this section of particular interest to the practical engineer, a number 


of typographical errors occur in identifying the graphs. 
Comparisons are made between theoretical calculations and experimental measurements obtained 
and the constructional difficulties involved in precise measuring systems are 


from various sources, 
discussed in considerable detail. In one reference in which this reviewer participated, however, it is 
noted that the author was incorrect in his statements describing the procedure. 

Chapter III is devoted to a general investigation of the mutual coupling between antennas in 


various geometric configurations. An analysis is made of some of the more common types of antennas 
such as coplanar arrays, parasitic elements, folded dipoles, V-antennas, asymmetrically driven antennas, 
etc 

Chapter IV is devoted to the general analysis of the essential properties of receiving and scattering 
antennas. The freespace patterns and gains of various types of linear radiators are taken up in chapter 
V. Chapter VI discusses the electromagnetic fields of various configurations of linear radiators in 
commonly-used arrays. 

Chapter VII is devoted to a study of the primary electromagnetic field and radiation characteristics 
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of antennas over a conducting earth. This portion could well have been omitted. The effect of the iono- 
sphere is not taken up in this book. 

In the last chapter, the antenna is formulated as a boundary-value problem with primary emphasis 
on the mathematical rather than the engineering approach. The radiation characteristics of hemi- 
spheroidal, conical, and cylindrical antennas are considered in detail. Recent works on large-angle conical 
antennas are not described in this book. 

The appendix contains tables of generalized sine and cosine integrals useful in the theoretical calcu- 
lations, a considerable number of problems based on the presented material, and a bibliography on 
antenna theory and measurements. 

In general, the author seems inclined to favor the more complex approaches to the problems. For 
example, a simpler method of analysis for the ground plane antenna described in the literature* has not 
been included. Taken as a whole, however, the book provides a good background and reference for both 
the mathematician and applied scientist interested in this particular field. The book will be of relatively 
less usefulness to the general antenna engineer working with more complicated feed systems and antennas 
of large cross section which are beyond the scope of the simplified arrangements described. 


O. M. Woopwarb, Jr. 


Hydrodynamics. By Hugh L. Dryden, Francis D. Murnaghan and H. Bateman. Dover 
Publications, Inc., New York, 1956. 634 pp. $2.50 (paper-covered). 


As is indicated on the reverse of the title page, this is a republication of Bulletin No. 84 of the 
National Research Council, which has long been out of print. 


Deformation and flow of solids. Edited by R. Grammel. Springer-Verlag, Berlin, Gottingen, 
Heidelberg, 1956. 


Work on plasticity of solids has developed in the past mainly along two lines: 

1. Physical metallurgists and solid state physicists have attempted to interpret the basic plastic 
properties of materials in terms of their atomistic structure; the creation of the dislocation theory 
represents the most impressive example of this type of activity. 

2. Mechanical engineers and applied mathematicians have attempted to establish certain g 
principles of a macroscopic nature for the deformation of solids; these principles permit, for example, 
computations of the plastic behavior of materials under more complicated cases of loading, if the behav- 


eneral 


ior under simple loading conditions is known. 

In both of these two branches of the theory of plasticity, much progress has been realized in the past 
ten years. However, there has been only little interaction between the two. Since it is felt by many 
that a close connection between the two branches would be fruitful for the development of either of 
them, the International Union of Theoretical and Applied Mechanics held a meeting in Madrid in 
1955 on the deformation of solids, with the express purpose of including all aspects of the subject 

The book contains the lectures given at this meeting and the contributions to the discussion. The 
main topics are: theory of dislocations, theory of fracture, mathematical theory of plasticity, non-linear 
elasticity, viscoelasticity and relaxation. Some of the articles are in the form of surveys and others in the 
form of representative original papers. In this way the book easily provides a certain understanding 
of the present status of even those parts of the theory of plasticity with which the reader may be less 
familiar. It therefore should be of value to anyone who is interested in the deformation of solids 

K. Litcx: 


Handbuch Der Laplace-Transformation. By Gustav Doetsch. Band III. Birkhauser 


Verlag, Basel and Stuttgart, 1956. 300 pp. $9.35. 
Volume 2 of this multi-volume work concluded with the treatment of ordinary differential equations 
and this third volume begins, appropriately, with a treatment of partial differential equations. After 


*“An Ultra-High-Frequency Antenna of Simple Construction,’”’ G. H. Brown and J. Epstein, Com- 


munications, July 1940, p. 3. 
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a brief chapter defining the boundary value and initial value problems, a detailed treatment of the ap- 
plication of the Laplace transform to second order hyperbolic, parabolic, and elliptic equations with 
constant coefficients is given. 

Briefly considered in the third chapter are the operational treatment of a variable coefficient equation 
in which (1) coefficients are independent of the transformation variable, and alternatively (2) the coeffi- 
cients may be linear in the transformation variable. There follow chapters dealing with compatibility 
and uniqueness and with Huygens’ and Eulers’ principles. 

Subsequent parts of the book are concerned with the application of the transform technique to 
difference equations and to integral equations whose kernel is necessarily of the convolution type, both 
for the finite and infinite domains. 

The last chapter deals with the finite Laplace transform. Finally, an historical survey and a detailed 
bibliography are provided. The book should be of considerable interest and value to those who want 
a careful detailed spelling out of the use of the Laplace transform in connection with the problems 


mentioned above. 
GeEorGE F. CARRIER 


Theoretische Hydromechanik. By N. J. Kotschin, I. A. Kibel, and N. W. Rose. Translated 
from the Russian into German by K. Krienes. Akademie-Verlag, Berlin. Volume 
I, 1954, XII + 508 pp.; Volume II, 1955, VIII + 569 pp. $11.50 each volume. 


The translation is of 1948 Russian edition. Volume I is concerned with the classical theory of the 
motion of an ideal fluid. The contents are as follows, the figures in brackets indicating the number of 
pages devoted to each topic; Chapter I, Kinematics of fluid motion: Deformation of a fluid drop (7): 
equation of continuity (15): kinematic character of irrotational and rotational flow (12). Chapter II, 
Basic equations of the hydrodynamics of an ideal fluid (36). Chapter III, Hydrostatics: Hydrostatic 
pressure (12): equilibrium of floating bodies (13). Chapter IV, the simplest motions of an ideal fluid: 
The integrals of Bernoulli and Cauchy (19): plane irrotational flow (14). Chapter, V, Vortex motion of 
an ideal fluid: The basic equations of vortex theory and the Helmholtz vorticity conservation theorems 

27): the determination of the velocity field from a given vortex and source field (31): the Kirmaén 
vortex street (30). Chapter VI, The plane problem for the motion of a body in an ideal fluid (107). 
Chapter VII, The three-dimensional problem for the motion of a body in an ideal fluid (41). Chapter 
VIII, Wave motions of an ideal fluid: Basic equations of wave theory (7): plane waves (67): three- 
dimensional waves (21): Jong waves (31). The apparently disproportionate length of Chapter VI is 
explained by the fact that it treats also of free streamlines, thin aerofoils, and (in some detail) the planing 
of a plate at the surface of water. The plate here considered is of infinite breadth, although the case of 
the plate of finite breadth was solved by A. E. Green as long ago as 1936. 

The subject is treated in great detail and there are some exercises which enhance the value of the 
work. On the other hand nothing like full use is made of vectors and the complex variable. Thus, for 
example, three pages are used to prove Kelvin’s theorem on the constancy of circulation and Cartesian 
presentation is used throughout the section on plane waves. It is, however, gratifying to see convincing 
expositions of some of the physical deductions concerning such diverse matters as generation of vorticity, 
trade winds, breezes, and ocean currents due to variable salinity, which flow from the theorem of Bjerknes. 
This theorem expresses the rate of change of circulation in a circuit in terms of the number of unit tubes 
defined by isobars and isobulks which thread the circuit. 

The contents of Volume II are as follows; Chapter I, Theoretical foundations of gasdynamics (by 
Kibel): The equations of gasdynamics (21): steady motion, the plane problem (143): steady motion, 
three-dimensional, including conical, flow (52): time dependent flow (24). Chapter II, The motion of 
viscous fluids (by Kotschin): Basic equations of motion of viscous fluids (46): exact solutions of the 
equations of motion (61): approximate solutions for small Reynolds number (42): approximate solutions 
for large Reynolds number, including boundary layer theory (86) and Oseen’s theory of vanishing 
viscosity (24). Chapter III, Elements of turbulence theory (by Kibel): Turbulence and instability (27): 
fully developed turbulence (17): mean values of the hydrodynamic quantities (9). 

In this volume also the treatment is detailed in some cases to the point of prolixity. Thus in Chapter 
II a good tensor notation could have reduced to about one the eleven pages which end up with the 
assumption that in a viscous fluid 3A + 24 = 0. It should be noted, and this applies to both volumes, 
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that 8 years have elapsed since the publication of the work in Russian so that it may be inferred that 
the treatise describes the state of the subject as it existed at least 10 years ago. In the case of Volume 
II therefore the reader must not be disappointed to miss the great advances which have been made in 
recent years. With this cautionary remark, it may be stated that the book gives a good and readable 
account of the classical parts of the subject matter listed above. 

L. M. Mrtnge-THomson 


Spheroidal wave functions. By J. A. Stratton, P. M. Morse, L. J. Chu, J. D. C. Little, 
and F. J. Corbato. The Technology Press of M.I.T., Cambridge, John Wiley & 
Sons, Inc., New York, and Chapman «& Hall, Ltd., London, 1956. xiii + 613 pp. 


$12.50. 


In the foreword to these tables of coefficients of expansions of spheroidal wave functions, P. M. Morse 
points out that while it is known that there are eleven systems of coordinates which permit separation 
of the scalar wave equation, only three can be conveniently used in practice since only in the cases of rec- 
tangular, circular cylindrical, and spherical coordinates have the corresponding functions been tabulated. 
The solutions of the wave equation for elliptic cylindrical coordinates lead to Mathieu functions which 
have been treated in some detail in the literature, and for which some tables are available. The book 
under review presents corresponding material for spheroidal wave functions. 

A section by Chu and Stratton discusses various expressions for these functions and mathematical 
relations between them. This is a reproduction of an earlier paper from the Journal of Mathematics and 
Physics (XX, 1941). The application of this material to the determination of the coefficients tabulated 
is presented in a second section by Little and Corbato. They discuss the numerical methods utilized, and 
emphasize the influence of the application of a high speed digital computer in carrying out this work. The 
computer was programmed and the output arranged to print the material in the form in which it appears 
in the book. This system avoids errors of computation and errors due to transcriptions between the 
calculated result and the printed tables. 

The tables give the coefficients of the expansions of both oblate and prolate spheroidal radial func- 
tions in series of spherical Bessel functions, and of the corresponding angular functions in series of 
associated Legendre functions. References to auxiliary tables of spherical Bessel and associated Legendre 
functions are given. Tables of the spheroidal wave functions themselves were not presented at this 
time due to the volume of material involved with the many subsidiary parameters. The coefficients 
tabulated permit numerical values of the functions to be obtained quite quickly. 


E. H. LEE 


Engineering analysis. A survey of numerical procedures. By Stephen H. Crandall 
McGraw-Hill Book Co., Inc., New York, Toronto, London, 1956. x + 417 pp. $9.50. 


This well-written volume, latest member of the Engineering Societies Monographs series, is concerned 
with the analysis of complex engineering problems and the methods appropriate to their numerical 
solution. It is addressed to engineers or engineering mathematicians at the graduate level, and is based 
on lecture courses given by the author at the Imperial College of Science and Technology and the 
Massachusetts Institute of Technology. By far the largest part of the book is devoted to a survey of 
important numerical techniques; however, the prior task of reducing a problem to mathematical form is 
also recognized and discussed. 

Perhaps the outstanding feature of the book is its organization, not according to a classification of 
numerical methods, but according to a classification of the problems to which they may be applied. 
Thus a given computational scheme, for example, relaxation, may be treated in several places. Three 
classes of problems are considered: equilibrium, eigenvalue, and propagation problems. The author 
devotes two chapters to each class, one dealing with systems with a finite number of degrees of freedom, 
the other with continuous systems. Each chapter contains roughly the following: (i) a derivation of 


several typical problems; (ii) a review of the salient mathematical characteristics of the type of problem 
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in question; (iii) a discussion of the most suitable techniques for obtaining numerical solutions, usually 
accompanied by examples worked out wholly or in part. Remarks concerning the convenience of particular 
methods for hand or machine computation are included, as are brief discussions of the possible errors 
involved. 

In the opinion of the reviewer this book is a valuable contribution to the literature, particularly 
since it is written from the point of view of the engineer or applied mathematician who seeks to make 
use of numerical computations, rather than from that of the numerical analyst, who may be chiefly 
concerned with the refinement of the methods themselves. 

Since this work is designed as a survey and not as an exhaustive treatise, the treatment may at times 
be somewhat condensed for the student entirely unfamiliar with the subject. Copious references to the 
existing literature, however, will guide the serious reader to more detailed accounts of those topics about 
which he requires more information. The book will be valuable as a reference to anyone concerned with 
the more practical aspects of numerical analysis, and should also be quite suitable as a text for a senior 
or graduate course in the field. Its worth in the latter connection is enhanced by the unusually large 
collection of exercises, some of which are drill problems, while others extend the ideas of the text. Most 
are accompanied by answers or hints as to their solution. The subject matter is well-indexed, and refer- 
ences to other works are provided by numerous footnotes, a selective bibliography, and an index of the 
authors cited. The physical appearance of the book is first rate, adhering to the standard set by previous 
members of the series. There are a number of minor errors, but the author has prepared a list of those 
which have come to his attention. 

WituraM E. Boyce 


Symposium on Monte Carlo methods. Edited by Herbert A. Meyer. Held at the University 
of Florida. Conducted by the Statistical Laboratory. Sponsored by Wright Air 
Development Center of the Air Research and Development Command. March 16 
and 17, 1954. John Wiley & Sons, Inc., New York, and Chapman & Hall, Ltd., 
London, 1956. xvi + 382 pp. $7.50. 


The proceedings of the symposium contains twenty papers covering a fairly wide range of topics 
including the generation of random numbers, general theory, and applications of Monte Carlo methods. 
Several papers are reviews written in a simple language for the novice whereas other are more specialized. 
Unfortunately, the papers are printed in the order of presentation at the symposium rather than in a 
logical order based upon content or difficulty. A reader unfamiliar with the subject could, however, 
understand nearly all the papers if he were to read them in a certain order starting with the reviews. 

The book also contains an excellent bibliography divided into three sections. Part I lists paper and 
abstracts on ‘Monte Carlo Proper’ including many rather inaccessible reports published by various 
laboratories. Part II lists references and abstracts on random digits and is meant to be fairly complete. 
Part III contains a selected list of articles and books (with some abstracts) mostly on related topics 
such as sampling. 

Despite the fact that there is no obvious attempt to organize the material as would be expected of a 
text book, the end result is actually better than many books supposedly written in an organized way. 
It is certainly a valuable and timely contribution to the literature. 


G. F. NEWELL 


Einfiihrung in die mathematische Statistik. By Dr. Leopold Schmetterer. Springer- 
Verlag, Wien, 1956. vi + 405 pp. $11.65. 


The author, noting the lack of any modern book on mathematical statistics in the German language, 
presents a fairly complete survey of developments of the last quarter century. Although labeled an 
introduction, it is not an elementary level book. The reader is assumed to be familiar with matrix theory, 
calculus and elementary set theory. With this as a starting point, the author concentrates on the subject 
at hand with a minimum of distraction. 

About a third of the book is an introduction to probability theory, primarily an outline of the modern 
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approach with many theorems stated without proof. This is followed by chapters on elementary sample 
theory, confidence regions, theory of parameter estimating, introduction to test theory, regression 
theory, introduction to nonparametric theory, and the classical Bayes method. There is also a German 
to English translation of most of the technical words. 

The style of writing is concise, easy to translate, and orderly. Although certain important topics may 
be left out, no topic is treated at excessive length. The book is well suited as a text book even for English 
speaking students preferably at the graduate level in American universities. 
G. NEWELL 


Elasticity, fracture and flow, with engineering and geological applications. By J. C. Jaeger. 
Methuen & Co., Ltd., London; John Wiley & Sons, New York, 1956. viii + 152 pp. 
$2.50. 


A thorough study of the mathematical theory of elasticity is often held to be the obvious first step 
in exploring the mechanics of deformable solids. In favor of this view speaks the fact that the math- 
ematical theory of elasticity sets standards of mathematical rigor that have not yet been, and probably 
never will be, attained in other branches of solid mechanics. On the other hand, prolonged exclusive ex- 
posure to the theory of elasticity may foster patterns of thinking that will make it unnecessarily difficult 
for the student to absorb the differing viewpoints of other branches of mechanics of solids. In recognition 
of this danger, several schools are currently experimenting with the idea of a general introductory 
course in mechanics of continua that is to be followed by more specialized courses on the various branches 
of mechanies of solids and fluids. The lack of suitable textbooks has severely handicapped these efforts. 
The present volume fills this gap for solids. 

The first chapter (48 pp.) is concerned with the analysis of stress, and finite and infinitesimal strain. 
In the second chapter (68 pp.), the author discusses the mechanical behavior of actual materials and 
introduces suitable mathematical models of elastic, viscous, plastic, and visco-plastic behavior. The 
third chapter (43 pp.) treats a number of typical elementary problems concerning equilibrium and 
motion of these solids. 

The exposition is clear and easy to follow. Throughout the book, the emphasis is on the basic assump- 


tions and the manner in which these affect the solutions. 
W. PRAGER 


Modern mathematics for the engineer. Edited by E. F. Beckenbach. McGraw-Hill Book 


Company, Inc., New York, 1956. xx + 514 pp. $7.50. 


The chapters of this book correspond to a series of invitation lectures given at the University of 
California, Los Angeles. The book consists of three parts. Part I is entitled Mathematical Models and 
contains the following chapters: Linear and Nonlinear Oscillations (S. Lefschetz), Equilibrium Analysis: 
The Stability Theory of Poincaré and Liapunov (R. Bellman), Exterior Ballistics (J. W. Green), Elements 
of Calculus of Variations (M. R. Hestenes), Hyperbolic Partial Differential Equations and Applications 

R. Courant), Boundary Value Problems in Elliptic Partial Differential Equations (M. M. Schiffer), 
lastostatic Boundary-Value Problems (I. S. Sokolnikoff). Part II is devoted to Probabilistic 
lems; it contains the following chapters: The Theory of Prediction (N. Wiener), The Theory of 
H. F. Bohnenblust), Applied Mathematics in Operations Research (G. W. King), The Theory 
amic Programming (R. Bellman), Monte Carlo Methods (G. W. Brown). Part III has the title 
putational Considerations and consists of the following chapters: Matrices in Engineering (L. A. 
Pipes), Functional Transformations for Engineering Design (J. L. Barnes), Conformal Mapping Methods 
(E. F. Beckenbach), Nonlinear Methods (C. B. Morrey, Jr.), What Are Relaxation Methods? (G. E. 
Forsythe), Methods of Steepest Descent (C. B. Tompkins), High-Speed Computing Devices and Their 





Applications (D. H. Lehmer). 
W. PRAGER 

























































apireonee 














